Short-term forecasting of wind speed and related electrical power

Short-term forecasting of wind speed and related electrical power

PII: S0038-092X(98)00032-2 Solar Energy Vol. 63, No. 1, pp. 61–68, 1998 © 1998 Elsevier Science Ltd All rights reserved. Printed in Great Britain 003...

143KB Sizes 4 Downloads 123 Views

PII: S0038-092X(98)00032-2

Solar Energy Vol. 63, No. 1, pp. 61–68, 1998 © 1998 Elsevier Science Ltd All rights reserved. Printed in Great Britain 0038-092X/98 $19.00+0.00

SHORT-TERM FORECASTING OF WIND SPEED AND RELATED ELECTRICAL POWER M. C. ALEXIADIS,* P. S. DOKOPOULOS,* H. S. SAHSAMANOGLOU** and I. M. MANOUSARIDIS*

* Power Systems Laboratory, Electrical and Computer Engineering Department, Aristotle University of Thessaloniki, 54006 Thessaloniki, Greece ** Meteorology and Climatology Department, Aristotle University of Thessaloniki, 54006 Thessaloniki, Greece Received 3 April 1997; revised version accepted 21 January 1998 Communicated by JOACHIM LUTHER

Abstract—Wind speed and the related electrical power of wind turbines are forecasted. The work is focused on the operation of power systems with integrated wind parks. Artificial neural networks models are proposed for forecasting average values of the following 10 min or 1 h. Input quantities for the prediction are wind speeds and their derivatives. Also, spatial correlation of wind speeds and its use for forecasting, are investigated. The methods are tested using data collected over seven years at six different sites on islands of the South and Central Aegean Sea in Greece. © 1998 Elsevier Science Ltd. All rights reserved.

than directly on power time series (Beyer et al., 1994) and this has also been adopted in the present paper. Time series of wind speed V are transformed into series of power P using manufacturers’ curves corrected for the air density. In the case of a wind farm the total power output can be reliably assessed by curves relating it with the wind speed at one or more reference points. To check the quality of the forecasting model the errors are compared with the ones given by the ‘‘persistent forecasting’’ of order p, denoted here as PF( p). There, the forecasted value T time ahead, V (t+T ), is in general equal to the f moving average of the last p values. Mostly, only the last value is used ( p=1). Earliest forecasting attempts referred to in the literature included a typical time series analysis that led to various Box–Jenkins ARMA ( p,q) models. The p most recent speeds and the q most recent forecasting errors were used as inputs. Whenever the sampling time t that was used for inputs, equalled the time ahead T, forecasting errors were similar or worse than the persistent (Fellows and Hill, 1990). The models seldom used the original time series as inputs (Contaxis and Kabouris, 1991), in most of the cases the data were first transformed in order to abort non-stationary or impose normal-type distribution (Hu et al., 1991; Daniel and Chen, 1991; Desrochers et al., 1986). Bossanyi (1985) and Beyer et al. (1994) used

1. INTRODUCTION

Wind turbines ( WTs) are frequently used as fuel savers especially in small autonomous electrical power systems with high operational costs. Therefore, an increased wind power penetration is of interest. However, an increase beyond a certain percentage (a few per cent) may considerably influence the power quality and the operation of power systems and may become hazardous for the conventional production means. Variations in wind power cause, in general, voltage and frequency fluctuations. Also, a sudden cut-off of wind power due to excessive wind speeds may cause unacceptable shocks in the conventional power units. Increased power penetration is possible if measures are taken concerning wind turbine control, unit commitment, reserve limits and dispatch of power units. This procedure may also affect the energy efficiency of the conventional power stations, since it affects the operation point of the power units. Prediction of wind power is therefore of importance for efficient load management and operation of the system. We assume here that electrical power output in a WT depends only on the amplitude of the local velocity vector of the air particles, which is function of time t, i.e. V(t) and is denoted here as ‘‘wind speed’’. It has been suggested that WT power forecast should be based on wind speed forecast rather 61

62

M. C. Alexiadis et al.

Kalman filters and artificial neural networks (ANNs), respectively, and suggested that the sampling time for the inputs should be shorter than the time ahead (e.g. t=10 s to 1 min for T=10 min). The errors reported were about 5–10% better than persistent. The above mentioned methods used only the historical data of the site under consideration, i.e. no values related to the physical phenomenon, as pressure profiles or relation of wind speed at different sites. This is probably the reason for any limitation in their efficiency, especially in long-range forecasting, e.g. one hour ahead. This has been recognised and alternative methods were proposed, based on relations of wind speeds among several sites (Chan et al., 1983; Schlueter et al., 1984, 1985). A forecasting formula was proposed on the basis that a disturbance in wind speed, travels from site A to site B, in a certain time and it is generally changed in its amplitude; its validity was reported to be limited to a certain wind speed profile (Schlueter et al., 1985). The present paper suggests models based on artificial intelligence and spatial relations of the wind speed so as to improve the efficiency of short-range forecasting. Investigations are supported by measurements at six sites on islands in the Aegean Sea, taken over the last seven years. 2. WIND CHARACTERISTICS

We realised that prior to forecasting it was useful to investigate the physical cause and the time series of the site under consideration. 2.1. Windy regime The local conditions may influence what we call the ‘‘wind profile’’, that is the speed and direction of the wind as a function of time. Wind velocity, as vector field, is determined by pressure differences (‘‘high’’–‘‘low’’) and the boundary conditions, i.e. the terrain. Pressure differences can also locally exist due to heat transfer between the earth’s surface and the air. The latter case is mostly valid near the sea and may induce the so-called ‘‘sea breeze’’. In our measurements, sea breeze was not significant. The area of the Central and South Aegean Sea is regarded as the most windy area in Greece. Prevailing winds are mostly NW to NE and are caused by the following synoptic situation that dominated during our measurements. During the cold months a ‘‘high’’ in the NW

region of the Balkan Peninsula is combined with a ‘‘low’’ in the East Mediterranean or in Cyprus. The synoptic situation during the warm sixmonth season differs from the previous in two ways. Firstly, the centre of low pressures is not dynamic, but thermal; secondly, the centre of this system is situated in the major region of Iraq and Saudi Arabia. The related winds are strong northern winds, called Etesian winds. 2.2. Measurements Instruments used were manufactured by NAG; a cup anemometer, a potentiometer driven wind direction indicator, a pressure and a temperature sensor. Values of 2 s were sampled on a PC and then processed through different filters to supply us with mean, medium, mean square or mean cubed values over periods of 1, 5 and 10 min or 1 h. 2.3. Pre-processing, investigation of time series The selection of the inputs is an important stage in building an efficient forecasting model. (1) Usually, a statistic analysis of the time series may show which previous values are mostly related with the future ones. The autocorrelation curve shows that the signal is nonstationary (Fig. 1). This was expected for the wind at the Aegean Sea, the average speed of which can move from 3 to 20 m s−1 in two days. Also there is a slight daily periodicity, usually in the summer ( Fig. 2). (2) The prevailing winds in the summer months, known as Etesians, are strong northern winds veering east during daytime. They show a slight decrease during the night (Fig. 2) and appear for 12–20 days every month. The meteorological services can forecast an Etesian wind, one day ahead. Etesians that can last five successive days are developed by pressure gradients of at least 1 hPa/100 km.

Fig. 1. Correlation of the next-hour average value with 60 previous values of various t samplings (1 min to 1 h). The curves are based on annual data for the island of Syros.

Short-term forecasting of wind speed and related electrical power

63

especially in complicated problems with nonsimilar inputs and many outputs. In our case, it was found that hidden layers did not improve the ANN’s performance, as was also suggested by Beyer et al. (1994). This was probably due to over-parameterisation of the model. 3.1. Ten-minute forecasting

Fig. 2. Average daily profile in Syros for the month of August. The solid line shows averages of 1 h and the dashed lines show the standard deviation. In all the islands studied, there wasn’t an impressive average daily profile. The influence of the clock time during the day seems mostly to be non-significant.

Therefore, the probable inputs of the model are as follows: (1) Recent wind speeds of that same site. (2) Reported speeds from neighbouring sites’ winds that are directed towards our site (spatial correlation). (3) The meteorological services’ global forecasts may be used as auxiliary inputs. They are usually given in the Beaufort scale with, relatively, great error. In this paper we used (1) and (2) as inputs while (3) as an input is used in a work to be reported in the future. 3. FORECASTING WITH ARTIFICIAL NEURAL NETWORKS

In this chapter forecasting results are reported for 10 min and 1 h ahead using as inputs only the historical data of the site under consideration. Both ARMA models and artificial neural networks (ANNs) were initially used. Best results were obtained using ANNs. ANNs have been a preferred technique in many forecasting or simulation applications for the last five years. An ANN consists of ‘‘neurones’’ structured in layers; an input layer, an output layer and between them one or more hidden (middle) layers. In each neurone the input applied is biased and then filtered, via a sigmoid-type function. In a typical fully connected ANN each neurone is linked with all the neurones of the previous layer and each link has its own ‘‘weight factor’’. The ANN parameters, i.e. weights and biases, are determined using the back-propagation method (Maren et al., 1990). The middle layer of the ANN usually plays a synthesis and filtering role for the inputs,

The average 10 min value shows a decreasing correlation with previous speeds ( Fig. 1). Therefore, it is reasonable to use the N most recent minute values as inputs. The minimum value of N, which is of interest (e.g. 5, 10, 15 or 30) has been determined after tests so as to give the least errors. As for as input processing is concerned, two methods were investigated. (1) the regular patterns method where speeds are used in their native form, i.e. as reported in m s−1; and (2) the ‘‘differenced’’ patterns method, suggested in the present paper, where the inputs V , i=1,…,N, and the output dif,i V of the patterns are expressed as dif,out follows: V =V −V , V =V −V dif,i i ave dif,out out ave

(1)

where N V = ∑ (V /N) ave i I=1 V is the ith input V is the output i out Comparing the two methods, we found that the second model appears superior in the following features: (1) errors are smaller in all cases, (2) the training set can be much smaller, (3) training takes a much shorter time (about 250 iterations instead of 2500). When regular patterns are used the training of ANN requires a wider data set. The model appears to be weak for cases that were not included in the training set or were badly assimilated due to their scarceness (e.g. strong winds of 25–30 m s−1). It is worth mentioning that the same model will have to satisfy a great range of values from 0 to 30 m s−1. For these reasons differenced patterns are preferred. They cope with the winds as deviations from the recent average value, grouping, in this way, all the cases ( low, medium and high winds) into one case. They did not show rogue errors, nor persistence or polarity of the errors.

64

M. C. Alexiadis et al.

In all cases referred in this paper the average forecasting error is defined as follows n ANN Error= ∑ |V −V |/n (2) i,real i,forecasted i=1 with V being the wind speeds in m s−1, and the i per cent improvement of the error is (Persistent Error−ANN Error) (Persistent Error)

×100% (3)

The same way, we define wind power error using P and P in kW. real forecasted Each time, two separate data sets were used: a training set to adjust the parameters of the ANN and a recall set to check the model’s efficiency. The recall set consisted of eight months (April to November of 1991). The same procedure was repeated for the three Cyclades Islands: Syros, Paros and Kea. For each island several N:1 models were tested, where N was the number of the most recent 1 min speeds used as inputs (N=5, 10, 15, 20, etc.). Table 1n shows that optimum N varies from site to site. Usually N=15 is enough. 3.2. Forecasting per hour The average wind speed of the next hour is related mostly to the 2–3 recent hours’ speeds (Fig. 1). Differenced and regular patterns as mentioned in the 10-min forecasting model, are also used here. It is further worth noticing that the persistence forecast is expected to be unsatisfactory, since variation of wind speed between two successive hours is large. Investigation of our models showed the following results: (1) Various time t averaged data were tested Table 1. Per cent improvement of the forecasted wind power error as compared to the persistent error for 10-min forecasting on Kea island (improvement of wind speed error in parentheses) Number of 1-min inputs

Regular patterns (%)

5 +7.033 10 +8.047 15 +8.076 30 +8.512 60 +8.645 120 +8.417 Optimum on Syros island 15 (15) 15 (15) +7.897 Optimum on Paros island 30 (15) 15 (15) +9.146

Table 2. Per cent improvement of average wind power error (pu) as compared to the persistent error for one-hour forecasting on Kea island (respective improvement of wind speed error in parenthesis) Number of taveraged inputs

Regular patterns (%)

Differenced patterns (%)

60 30 20 6 6 5

+12.87 +12.80 +12.76 +11.52 +5.529 −1.536

+13.70 (+12.39) +13.86 (+12.24) +13.89 (+12.88) +12.70 (+11.82) +6.697 (+6.917) −0.2389 (−0.352)

1 min 2 min 3 min 10 min 30 min 1h

(+9.212) (+9.528) (+9.493) (+8.211) (+3.316) (−2.101)

as inputs (e.g. 15, 10, 5, or even 1 min). Inputs out of short sampling time t are suggested because there are many of them, helping the robustness of the model and they are all highly correlated with the output ( Fig. 1). For instance 20 inputs of three-minute values or 30 inputs of twominute values may be used. (2) Differenced patterns give smaller errors than regular patterns. The improvement they bring is increased especially when the only input data disposed are of long sampling time (t≥5 min). Results are presented in Table 2. Characteristic curves of measured and forecasted wind speed and WT power time series are shown in Figs 3 and 4, respectively.

Fig. 3. Measured (solid line) and forecasted wind speed (m s−1), for one-hour averages. The model used is the optimum one of differenced patterns (shown in bold in Table 2).

Differenced patterns (%)

(+6.298) (+6.834) (+7.018) (+7.034) (+7.036) (+7.021)

+8.474 +9.321 +9.322 +9.536 +9.460 +8.948

(+6.171)

+8.082 (+6.808)

(+8.800)

+9.536 (+9.041)

(+8.114) (+9.158) (+9.134) (+9.008) (+8.696) (+8.141)

Fig. 4. Measured (solid line) and forecasted WT power (pu) for the data of Fig. 3.

Short-term forecasting of wind speed and related electrical power

65

4. SPATIAL CORRELATION MODELS

The application of spatial correlation models will be now examined in the area of the Central Aegean Sea (Cyclades Islands). The validity of the models will be tested using data from three sites where wind turbines were installed. The sites are at altitudes of 400, 547 and 724 m on the islands of Syros, Kea and Paros with nominal powers of 100, 60 and 100 kW, respectively. These three islands form a triangle with sides of about 60, 50 and 105 km (Fig. 5), distances far larger than the ones referred to in previous studies such as Schlueter et al. (1985). In the following, two methods will be investigated: the first one is based on the crosscorrelation curve between the speeds at the two sites (Schlueter et al., 1985); the second method, is suggested by the authors as a new spatial correlation predictor (SCP). 4.1. Forecasting using the cross-correlation curves (constant delay method) Forecasting models suggested in the literature are often based on cross-correlation curves of wind speeds at two sites A and B (Schlueter et al., 1984, 1985). Winds are blowing from A to B, i.e. there is a speed component directed from A to B. These models use as constant time delay the interval Dt at which cross-correlation becomes maximum and require that the peak correlation coefficient r(Dt) is high (≥=0.8). The following formula is then used to forecast future speed at B, V(B, t+Dt) from the present speed at A, V(A, t) using the average values m , m and standard deviations s , s A B A B V(B, t+Dt)=m +[r(Dt)s /s ][V(A, t)−m ] B B A A (4) Such a model was tested for the period of the summer months for the Syros (A) and Paros (B) islands. This case seemed to satisfy the requirements of the method since the crosscorrelation peak was high (0.85 to 0.9) and the same event (strong northern winds) was repeated on the majority of the days. Nevertheless errors were far worse than persistent. The main drawbacks and restrictions of the model are described as follows: Even when there is an enormous speed difference between the two islands a shape similarity in wind profiles is observed. So, for sudden speed changes, i.e. high speed derivatives occur at A, analogous changes appear later at B and

Fig. 5. Sites on Cyclades Islands, Greece, where the spatial correlation predictor was applied.

are, in general, successfully forecasted. Though, often, the model fails to bring the forecasted fluctuations on the right speed level. An error correction method would be indispensable to solve this problem. The constant time delay is often wrong; it appears that certain events arrive an hour earlier or later than expected. The model is bound to a certain case—the average case. Hence, for large deviations from that case the general picture of forecasting is disappointing and the model is not useful at all. 4.2. The proposed spatial correlation predictor (SCP) We consider two equally levelled points A and B close enough to each other. The points are at a sufficient height above sea level, in a layer where the stream lines are parallel and the vertical component of wind velocity vanishes. If the wind is blowing from A to B and at a time t a considerable change dV occurs in the velocity V(A, t), this change will appear at B after a time delay Dt. Let us denote the velocities at A and B and the propagation velocity of the event as V(A, t), V(B, t) and V (t), respectively. F Since stream lines are parallel we may assume that V(A ,t)=V(B, t+Dt)=V (t) (5) F Dt=(AB)V (t)V−2 (t)=(AB)V(A, t)V−2(A, t) F F (6) where (AB) is the vector of the distance (AB) and (AB) V(A, t) denotes the inner product, i.e. Dt is the projection of (AB) on V (t) divided by F the propagation speed. We can determine the direction of V (t) using the pressure measureF

66

M. C. Alexiadis et al.

ments at three points and referring them to sea level or any common level (Holton, 1979). Equation (5) implies that in a homogenous velocity field any changes of wind propagate with a speed equal to the wind velocity. Now let us consider any two points A∞ and B∞ below A and B, respectively. We would expect the velocities at A∞ and B∞ to correspond to the speeds at A and B as they are modified locally by the terrain, the altitude and thermal effects. In these cases linear relations were found adequate to estimate Dt and forecast wind speed at B∞ from its value at A∞ as follows: V (t)=V(A, t)=a +a V(A∞, t) (7) F 0 1 V(B∞, t+Dt)=b +b V(A∞, t) (8) 0 1 where a , a , b and b are constants to be 0 1 0 1 determined using a method described in the next chapter, Dt is the translation time which corresponds to the speed V (t) according to F Eqs. (6) and (7). 4.3. Implementation of SCP The model was tested using data from Syros and Paros, as points A and B. Two separate data sets, both from summer months, were used: one to optimise the model and one to check it. The following methods are advised for the parameters determination. Coefficients b and b are used to adjust wind 0 1 the speeds of A∞ to the wind characteristics of B∞ (see also eqn (4)), that is they are selected to avoid permanent under- or over-prediction. Their values should not be considered proper until the model’s average error becomes unbiased. Coefficients a and a specify the relation of 0 1 local and upper winds at site A. A false estimation of them would lead the model to forecast the oncoming events earlier or later than it should. Such a mistake can be easily traced by checking the forecasting graphs of a few characteristic cases. The wind speed sequences may, in general, show the following three states: (1) speed may flutter widely around a rather steady average level; (2) a sudden change may occur, i.e. a vast rise or fall of the average level; (3) a steady trend may appear for a long time interval. The SCP model is mostly appropriate for cases (2) and (3). These cases are not only the most impressive ones but also the most characteristic to judge the model’s efficiency.

Nevertheless, case (1) which covers a great percentage of time may lead to successive errors of the same polarity. This problem is eliminated by adding an error correction term to the forecast, i.e. adding a quantity c · e(t−1), where c is a parameter (here estimated as −0.75), and e(i−1) is the latest error of the non-corrected model. When the model is used for forecasting the wind speed at sites where no measurements are taken, an error correction term can not be used. The resulted time delays varied from 10 min to 5 h (Fig. 6) depending on the wind speed and, especially, the wind direction. Using this model for 10-min forecasting was not practical for such long distances and long delay times. Forecasting for one or more hours ahead was made using 10-min averaged speeds from Syros and translating them to Paros according to the model. The translated values were then averaged to give one-hour averages. Since, every 10 min a different time delay is calculated it appears that (1) There may be gaps in the forecasted time series. (2) Latest fast winds may overtake slow winds, that is, they may substitute a previous forecast for time t+Dt that was based on low winds. (3) Not always the same sort of forecast is possible. At present time t, you may have a forecast for the next 1, 2 or 3 h, or no forecast at all. So, an actual comparison with the other formerly-mentioned models of this paper is not easy for all cases. Nevertheless, in practice these problems are not that important. Usually, there is a long interval of the same event, e.g. 8 h, during which there is a varying speed magnitude and a steady

Fig. 6. Frequency of time delay Dt calculated by the spatial correlation predictor.

Short-term forecasting of wind speed and related electrical power

67

Table 3. Per cent improvement of average wind power error (pu) of the SCP model as compared to the error of the constant delay method for one-hour forecasting on Syros island (respective improvement of wind speed error in parenthesis) a 0

a 1

b 0

0 0 0.80 0.80 0.80

1 1 1 1 1

0 0 1.10 1.10 1.10

b 1

c

Improvement (%)

1 1 1.15 1.15 1.15

0 −1 0 −1 −0.75

−74.33 (−12.93) −9.526 (19.92) −6.192 (9.238) 13.42 (30.67) 16.08 (32.14)

or smoothly changing direction. So, at each time you normally dispose e.g. 5 to 8 10-min forecasts for the next 90 min. For a reference site, e.g. Syros in this case, the model can only be used for a limited duration of the year (about 70% of the time), this depends on the wind direction which is known. Though, it should be noted that the remaining 30% mostly consists of low, non-event, noncorrelated wind speeds. Table 3 presents cases of next-hour forecasting only, in order to compare it with the constant delay method (in which Dt was 50 min), and also shows the sensitivity of the model over its parameters. An example of SCP implementation using the optimum model of Table 3, is shown in Figs 7 and 8. The winds blowing on that day were NW to N and the SCP suggested delay times from 1 to 4 h according to the speed magnitude. 5. CONCLUSIONS

Forecasting of wind power may be satisfactory by using artificial neural networks (ANNs). The differences of wind speeds from their moving average should be preferred as inputs.

Fig. 7. Forecast of wind speed using the spatial correlation predictor. (1) Measured speed at B; (2) measured speed at A; (3) forecasted speed at B, SCP method; (4) forecasted speed at B, constant delay method.

Fig. 8. Forecast of WT power using the spatial correlation predictor for the data of Fig. 7. (1) Measured power at B; (2) forecasted power at B, SCP method; (3) forecasted power at B, constant delay method.

The sampling time of the input values should be shorter than the time ahead of forecasting. Also, a model of spatial correlation is proposed for predicting wind speed at one site based on measurements at another site. Its behaviour has been tested with satisfactory verification using data collected over seven years. Acknowledgements—This project was financially supported by the General Secretariat of Research and Technology of the Greek Ministry of Development.

REFERENCES Beyer H. G., Degner T., Hausmann J., Hoffmann M. and Rujan P. (1994) Short term prediction of wind speed and power output of a wind turbine with neural networks. EWEC ’94, 5th European Wind Energy Association Conference And Exhibition, Conference Proceedings, Oral Sessions, Vol. 1, Thessaloniki, Greece, pp. 349–352. Bossanyi E. (1985) Short term stochastic wind prediction and possible control applications. In Proceedings of the Delphi Workshop on ‘‘Wind Energy Applications’’, Greece, pp. 66–79. Chan S. M., Powell D. C., Yoshimura M. and Curtice D. H. (1983) Operations requirements of utilities with wind power generation. IEEE Transactions on Power Apparatus and Systems 102, 9, 2850–2860. Contaxis G. and Kabouris J. (1991) Short term scheduling in a wind/diesel autonomous energy system. IEEE Transactions on Power Systems 6, 3, 1161–1167. Daniel A. R. and Chen A. A. (1991) Stochastic simulation and forecasting of hourly average wind speed sequences in Jamaica. Solar Energy 46, 1, 1–11. Desrochers G., Blanchard M. and Sud S. (1986) A MonteCarlo simulation method for the economic assessment of the contribution of wind energy to power systems. IEEE Transactions on Energy Conversion 1, 4, 50–56. Fellows A. and Hill D. (1990) Wind and load forecasting for integration of wind power into a meso-scale electrical grid. European Community Wind Energy Conference, Proceedings of an International Conference, Madrid, Spain, pp. 636–640. Holton J. R. (1979) An Introduction to Dynamic Meteorology. Academic Press, New York, p. 58. Hu J., Beans J. and Auslander D. M. (1991) Start/stop control of fixed-pitch wind energy turbines. Solar Energy 46, 1, 29–40. Maren A., Harston C. and Pap R. (1990) Handbook Of

68

M. C. Alexiadis et al.

Neural Computing Applications. Academic Press, San Diego, California, pp. 85–105. Schlueter R. A., Park G. L., Bouwmeester R., Shu L., Lotfalian M., Rastgoufard P. and Shayanfar A. (1984) Simulation and assessment of wind array power variations based on simultaneous wind speed measurements.

IEEE Transactions on Power Apparatus and Systems 103, 5, 1008–1016. Schlueter R. A., Sigari G. and Costi A. (1985) Wind array power prediction for improved operating economics and reliability. IEEE Transactions on Power Apparatus and Systems 104, 7, 137–142.