SW—Soil and Water: Application of the Simulated Annealing Method to Agricultural Water Resource Management

doi:10.1006/jaer.2001.0723, available online at http://www.idealibrary.com on SW*Soil and Water

Application of the Simulated Annealing Method to Agricultural Water Resource

Management

Sheng-Feng Kuo; Chen-Wuing Liu; Gary P. Merkley

Department of Civil Engineering, National Ilan Institute Technology, Taipei, Taiwan 114, ROC; e-mail:kuosf@chiseng.org.tw Department of Agricultural Engineering, National Taiwan University, Taipei, Taiwan 106, ROC; e-mail of corresponding author:

lcw@gwater.agec.ntu.edu.tw

Department of Biological and Irrigation Engineering, Utah State University, Logan, UT 84341, USA; e-mail:merkley@cc.usu.edu

(Received 5 April 2000; accepted in revised form 16 March 2001; published online 17 July 2001)

This work presents a model based on the on-farm irrigation scheduling and the simulated annealing (SA) optimization method for agricultural water resource management. The proposed model is applied to an irrigation project located in Delta, Utah of 394)6 ha area for optimizing economic pro"ts, simulating the water demand and crop yields and estimating the related crop area percentages with speci"ed water supply and planted area constraints.

The application of SA to irrigated project planning in this study can be divided into nine steps: (1) to receive the output from the on-farm irrigation scheduling module; (2) to enter three simulated annealing parameters; (3) to de"ne the design &chromosome' representing the problem; (4) to generate the random initial design &chromosome'; (5) to decode the design &chromosome' into a real number; (6) to apply constraints; (7) to apply an objective function and a "tness value; (8) to implement the annealing schedule by the Boltzmann probability; and (9) to set the &cooling rate' and criterion for termination.

The irrigation water requirements from the on-farm irrigation scheduling module are: (1) 1067)9, 441)7, and 471)8 mm for alfalfa, barley and maize, respectively, in one unit command area; and (2) 1039)5, 531)4, 490)9, and 539)4 mm for alfalfa, barley, maize and wheat, respectively, in the other unit command area. The simulation results demonstrate that the most appropriate parameters of SA for this study are as follows: (1) initial simulation &temperature' of 1000; (2) number of moves equal to 90; and (3) &cooling rate' of 0)95.

 2001 Silsoe Research Institute

1. Introduction

Agricultural water resource planning or irrigation planning can be described as a process for simulating complex climate}soil}plant relationships, applying mathematical optimization techniques to determine the most bene"cial crop patterns and water allocations. Such a determination can be non-trivial when large irrigated areas with signi"cant crop diversi"cation are considered, especially with the typical temporal and volumetric re-strictions on water supply. A computer-based model to simulate the climate}soil}plant systems with a new mathematical optimization technique could be an e!ective tool to help irrigation planners to make sound

decisions prior to each crop season. Maidment and Hutchinson (1983) stated that irrigation water manage-ment models may be classi"ed into two types: (1) demand simulation models, and (2) economic optimization mod-els. Demand simulation models pertain to the cli-mate}soil}plant system, and can be used to deduce the amount and timing of irrigation needed to ensure ad-equate crop growth. Economic optimization studies re-late the cost of irrigation to the bene"ts derived from increased crop productivity, among other possible fac-tors, to determine the economically optimal patterns of crops and irrigation water application.

Irrigation scheduling is a basic component of agricul-tural water resource planning. Many existing models

Notation a empirical coe$cient

aG, bG minimum and maximum values of de-_{coded decimal} A crop planted area, ha

AH% crop planted area within command

area, %

AHF? crop planed area within command area,_ha ASA? area of each unit command area, ha A  minimum percentage area of crop with-_{in command area, %} A  maximum percentage area of crop

with-in command area, %

dG depth of irrigation water, mm dL maximum net depletable depth, mm_{D soil moisture depletion, mm} D? soil allowable depletion, mm

D  soil maximum allowable depletion, mm_{E energy during annealing scheduling,} di-mensionless

EA, E? conveyance and application_{coe$cient, %}

EKMTC project bene"t at current move during annealing scheduling, $

E2 evapotranspiration, mm

E2M daily reference crop evapotranspiration,_{mm day\} E2A potential crop evapotranspiration,

mm day\

E2A? actual crop evapotranspiration, mm day\

daily reference crop evapotranspiration at each stage, mm day\

E2 AQR?EC potential crop evapotranspiration at each stage, mm day\

E2 A? QR?EC actual crop evapotranspiration at each_{stage, mm day\} DE change of project bene"t from current

and previous moves, $

fQC?QML cumulative seasonal in"ltration, mm FGH fertilizer cost of the jth crop in the ith_{command area, $ ha\}

i, j command area and crop index k decision variable

K? soil moisture stress coe$cient KA@ basal crop coe$cient

KRA@ basal crop coe$cient at day t

KQ coe$cient for evaporation rate from_{a wet soil surface} KW crop yield response factor

KWQR?EC crop yield response factor at current_{growth stage} LGH labour cost of the jth crop in the ith_{command area, $ ha\}

mG substring length M number of moves

M available soil moisture, mmm\ N number of command areas within

ir-rigated project

NA number of crops within command area OGH operation cost of the jth crop in the ith_{command area, $ ha\} PGH unit price of the jth crop in the ith

command area, $ ha\ PP Boltzmann probability

QBCK cumulative crop water demand in_{command area, m} QQSN available water supply for command area,_m

QGH cumulative water requirement of the jth_{crop in the ith command area, m} r random number

R? extraterrestrial radiation, mmday\ RX root depth, mm

s summation identi"er for substring length SGH seed cost per hectare of the jth crop in_{the ith command area, $ ha\}

t day of year

t, tL Julian days at the beginning and end of_{the crop growth stage} tB time required for soil surface to dry_{after irrigation or rainfall, days} tU time in days since wetting due to irriga-_{tion or rainfall, days} ¹ daily air temperature,C

¹ , ¹  maximum and minimum daily temper-_{atures, 3C} ¹Q? simulation &temperature' during cooling_{schedule, dimensionless} ¹LCU, ¹MJB simulation &temperatures' at the end_{and beginning, dimensionless}

= _{unit price of irrigation water, $ m}\ x decoded decimal

>?K crop yield reduction

>?QC?QML relative crop yield reduction due to_{in"ltration over the entire season} >?KQC?QML relative crop yield reduction due to_{water stress over the entire season} >?KQR?EC relative crop yield reduction due to_{water stress at each stage} >GH yields per hectare of the jth crop in the_{ith command area, ton ha\}

a &cooling' rate

hR soil moisture at the tth day

determine on-"eld water demands based on climate}soil}plant systems. Hill et al. (1982) developed the crop yield and soil management simulation model (CRPSM) to estimate crop yield as a function of soil moisture content, crop phenology and climate during the growing periods. Keller (1987) developed the unit com-mand area (UCA) model based partly on the concepts of the CRPSM model. The UCA model consists of two modules: the on-"eld module for water allocation and distribution; and the "eld and weather generation mod-ule. Prajamwong (1994) developed the command area decision support model (CADSM) with three main sub -models: (1) weather and "eld generation; (2) on-"eld crop}soil water balance simulation; and (3) water alloca-tion and distribualloca-tion. Smith (1991) developed the CROP-WAT computer program to calculate crop water require-ments and irrigation requirerequire-ments from climatic and crop data.

Simulated annealing is a stochastic computational technique derived from statistical mechanics for "nding near globally solutions to large optimization problems (Davis, 1991). The mathematical theory behind simulated annealing can be explained by the theory of Markov chains (Aarts et al., 1985; Otten and Ginneken, 1989) and in#uenced by the following three operators: (1) initial simulation &temperature'; (2) the number of moves to allowable rearrangements of the atoms within each temperature; and (3) the &cooling rate' to decrease the &temperature'. In a mathematical context, these three operators are the required parameters in the simulated annealing method. Kirkpatrick et al. (1983) were the "rst to propose and demonstrate the application of simula-tion techniques from statistical physics of combinasimula-tional optimizations. The mathematical theory to perform the idea of simulated annealing can be obtained using the theory of Markov chains (Laarhoven & Aarts, 1987). Bohachevsky et al. (1986) stated that the advantage of the simulated annealing method is the ability to migrate through a sequence of local extremes in search of the global solution and to recognize when the global extremum has been located.

It is interesting to review some papers in which simulated annealing has been applied to water resource management and irrigation scheduling (Dougherty & Marryott, 1991; Marryott et al., 1993). Dougherty & Marryott (1991) applied simulated annealing to three problems of optimal groundwater management: (1) a dewatering problem; (3) a contamination problem; and (4) contaminant removal with a slurry wall. Furthermore, they stated that "ve elements are needed to apply simulated annealing to a particular optimization prob-lem: (1) a concise representation of the con"guration of the system decision variables; (2) a scalar cost function; (3) a procedure for generating rearrangements of system;

(4) an annealing schedule; and (5) a criterion for ter-minating the algorithm. It was concluded that simulated annealing had the potential for solving groundwater management problems and that because the application of simulated annealing to water resources problems was new and its development is immature further perfor-mance improvements could be expected. Walker (1992) applied the simulated annealing method to a peanut growth model for optimization of irrigation scheduling. The peanut growth model was "rst applied to determine the days to irrigate and the amount of irrigation during the season. Later, simulated annealing was implemented in the peanut model. The general procedures of this study can be summarized as follows: (1) an initial vector with a "xed 10-day irrigation schedule was chosen from plant-ing to harvest to begin simulation for each year from 1974 to 1991; (2) the peanut model was run and a gross yield was obtained based on the initial vector selected; (3) a new vector of days was generated by the random number generator and a simulated crop yield was cal-culated; and (4) the Boltzmann distribution probability with generated random number was used to make the decision whether to update the irrigation days or not.

2. Model development

This study focuses mainly on developing an irrigation and planning model to simulate an on-farm irrigation system, and optimize the allocation of the irrigated area to alternative crops for maximum net bene"t by the use of a customized simulated annealing method. Therefore, this work develops a model based on the on-farm irrigation scheduling and simulated annealing method to support the agricultural water resource planning and management. The model consists mainly of six basic modules: (1) a main module to direct the running of the model with pull-down menu ability; (2) a data module to enter the required data by a user-friendly interface; (3) a weather generation module to generate the daily weather data; (4) an on-farm irrigation scheduling mod-ule to simulate the daily water requirement and relative crop yield; (5) a simulated annealing module to optimize the project maximum bene"t; and (6) a results module to present results by tables, graphs and printouts.

Six basic data types are required for the model: (1) project site and operation data; (2) command area data; (3) seasonal water supply data; (4) monthly weather data; (5) soil properties data; and (6) crop phenology and economic data. Herein, the weather generation module is adopted from CADSM (Prajamwong, 1994) to generate daily reference crop evapotranspiration and rainfall data based on the monthly mean and standard deviations data. The on-farm irrigation scheduling module receives

Fig. 1. The framework and logic employed in the agricultural water resource decision support model; E2, evapotranspiration the basic project data and generated daily weather data

to simulate the on-farm water balance. The daily simula-tion procedure includes three programming loops: (1) number of command areas within the simulated irriga-tion project, (2) number of crops within each command area, and (3) number of days from planting to harvest for each crop type. The daily simulation begins from the "rst command area in the project, the "rst crop within the command area, and the "rst Julian day for each crop. The procedure continues until all crops in each com-mand area and all comcom-mand areas in the irrigation project are processed. The output from this module in-cludes relative crop yield and crop irrigation water re-quirements. Both outputs are the required inputs for the following simulated annealing optimization module. Figure 1 presents the framework and logic employed in the irrigation decision support model.

3. On-farm irrigation scheduling

These on-farm irrigation scheduling processes deal with the daily water balance to estimate relative crop yield and irrigation water requirements. The Julian day of planting for each crop type is calculated based on speci"ed crop planting dates. Therefore, the Julian day at harvest is the sum of a crop planting day and cumulative

days from each of the growth stages. The daily simulation begins from the "rst command area in the project, the "rst crop within the command area, and the initial Julian day for each crop. This procedure continues until all crops in each command area and all command areas in the irrigation project are processed. The relative crop yield and irrigation water requirements are the return values from this module. The results are sub-sequently sent to the simulated annealing optimization module.

The method developed by Keller (1987) and Prajam-wong (1994) was used in this study to generate the daily weather data based on the mean monthly and standard deviation data. A normal distribution is assumed for generating the daily crop reference evapotranspiration and air temperature. A log-normal distribution is as-sumed for generating daily precipitation. Two important control factors are necessary for generating the daily weather data: (1) the arid probability; and (2) the random sowing date. The arid probability controls the aridity of the year and the random sowing date a!ects the sequence of the generated data. Based on the generated daily weather data, the Hargreaves equation (Hargreaves et al., 1985) was used to calculate the reference crop evapotran-spiration:

where: E2M denotes the (grass) reference crop evapotran-spiration in mm day\; R? represents the extraterrestrial radiation in mm day\; ¹ is the mean daily air ature in 3C; ¹  denotes the maximum daily air temper-ature in 3C; and ¹  represents the minimum daily air temperature in 3C.

The basal crop coe$cient KA@ represents the e!ects of the crop canopy on evapotranspiration and varies with time of year. The rate of change of the basal crop coe$-cient with time can be approximated as a linear increase (or decrease), as expressed in the following equation (Prajamwong, 1994):

KRA@"KQR?EC\

A@ #

(t!tQR?EC\);K_{tQR?EC!tQR?EC\}QR?ECA@ !KQR?EC\A@ (2a) and,

tQR?EC\)t)tQR?EC (2b)

where: KRA@ denotes the basal crop coe$cient for day t; KQR?EC

A@ represents the basal crop coe$cient at the current stage; tQR?EC\ is the "rst day of current crop stage; tQR?EC denotes the "rst day of the next crop growing stage; and t is the day of year.

The daily reference crop evapotranspiration E2M is used to calculate the potential E2A and actual crop evapotranspiration E2A?, as given in Eqns (3) and (4), respectively,

E2A"(KA@#KQ)E2M (3)

E2A?"(KA@K?#KQ)E2M (4) The soil moisture stress coe$cient K? and the coe$-cient for evaporation rate from a wet soil surface after irrigation or rainfall KQ are given by

K?"ln [100(hR!hUN)/(hDA!hUN)#1]ln(101) (5) KQ"(1!KA)





tB





U is the time in days since wetting due to irrigation and/or rainfall; and tB de-notes the time in days required for the soil surface to dry after an irrigation and/or rainfall event.

For on-demand irrigation scheduling, irrigation should be performed when the soil moisture depletion D initially exceeds the allowable depletion D?. The required amount, or application depth, dG in mm for a given irrigation, and allowable depletion D? in mm, can

be mathematically described by Eqns (7) and (8), respectively,

dG"_EAE?D (7)

D?"(hDA!hUN)RXD  (8) where: D represents the soil moisture depletion in mm; EA is the conveyance coe$cient; E? denotes the water application e$ciency; RX represents the root depth of the crop in mm; and D  is the maximum allowable soil water depletion in mm.

For each crop type, the cumulative water requirement in a growing season is the sum of the irrigation applica-tion depths at each time during the growing season. The cumulative water requirement for each command area is the sum of seasonal crop irrigation water requirements within the command area. Finally, the cumulative irriga-tion water requirement for the project is the sum of the water requirements of each command area within the project.

The amounts of in"ltration and runo! are calculated based on the irrigation water or e!ective rainfall multi-plied by the percentage of deep percolation and runo! due to irrigation and rainfall. The model user enters percentage values of in"ltration and runo!. The cumulat-ive amount of in"ltration is used to calculate the crop yield reduction due to waterlogging.

Two factors in#uence the relative crop yield: (1) the water stress due to insu$cient water for crop evapotran-spiration; and (2) waterlogging due to in"ltration, pro-duced by over-irrigation and or precipitation. Although the percentage of relative crop yield starts at 100% at the beginning of a growing season, the value can be reduced to less than 100% if there is any water stress or waterlog-ging during the growing season.

The relative yield reduction due to water stress is calculated at the end of each growth stage based on the ratio of cumulative potential crop evapotranspiration E2MQR?EC, and actual crop evapotranspiration E2AQR?EC in each stage. The relationships can be described by the following equations (Prajamwong, 1994):

>?KQR?EC"1!KWQR?EC





actual crop evapotranspiration at the end of the stage; E2AQR?EC denotes the potential crop evapotranspiration at the end of the stage; t and tL represent the Julian days at the beginning and end of the stage; and E2A? and E2A are daily crop potential and actual evapotranspiration in mm day\, respectively.

The minimum value of >?KQR?EC at each growth stage was chosen to be representative of the relative yield reduction due to water stress over the entire season >?KQC?QML as given by

>?K QC?QML" Min(>?K ; >?K ;2, ; >?KQR?EC) (12) The cumulative in"ltration within the root zone will reduce soil aeration due to waterlogging and in#uence the crop yield. Based on the only consideration of total in"ltration during the crop growth period, the relative yield reduction due to waterlogging is calculated at the end of the season based on the ratio of cumulative total in"ltration fQC?QML and the maximum net depletable depth dL in the root zone. These relationships can be represented by the following equations (Prajamwong, 1994):

>?QC?QML"1!a



dL"D MRX (14)

where: >?QC?QML denotes the relative yield reduction due to in"ltration over the entire season; a is the empirical coe$cient; D  represents the maximum allowable de-pletion (fraction); M is the available soil moisture in mm m\; and RX denotes the maximum root depth in m. The product of relative yield reduction due to water stress over the entire season >?KQC?QML and relative yield reduction due to waterlogging over the entire season >?QC?QML is the "nal value of relative crop yield at the end of the growing season.

4. Implementation of simulated annealing 4.1. Simulated annealing model

Simulated annealing (SA) has recently been applied to functional optimization problems. Functional optimiza-tion problems can be described as &real-world' problems with an objective of obtaining the minimum or maximum global values within speci"ed constraints. For decision support in irrigation project planning, this &real-world' problem attempts to obtain the optimal crop area-allo-cated values to maximize the bene"t of an irrigation project, given various constraints (e.g. maximum and minimum planted areas by crop type and maximum volume of water supply). The SA module has been imple-mented with the on-farm irrigation scheduling module

to maximize the project bene"t. The computational procedure of the SA module can be divided into the following steps: (1) to receive the output from the irriga-tion scheduling module; (2) to enter simulated annealing parameters through an user interface; (3) to de"ne the design &chromosome' to represent the problem; (4) to generate the random initial design &chromosome'; (5) to decode the design &chromosome' into a real number; (6) to apply constraints; (7) to apply an objective function and a "tness value; (8) to implement the annealing sched-ule by the Boltzmann probability; and (9) to set the &cooling rate' and criterion for termination. Figure 2 shows the #owchart of simulated annealing module. The following sections provide descriptive details about each of the steps.

4.2. Data requirements

Three parameters must be speci"ed for the simulated annealing module as follows.

(1) Initial simulation &temperature'

Following the steel industry analogy, the initial simu-lation &temperature' refers to the initial temperature for &annealing' in the model. The simulation &temper-ature' is gradually decreased depending on the simu-lation &cooling rate'. Also, the initial simusimu-lation &tem-perature' will in#uence the Boltzmann probability that dominates the annealing schedule.

(2) Number of moves

The number of moves is the allowable time for re-arrangement of the atoms to a lower energy state within each temperature value. Certainly, a higher number of moves will have a higher opportunity to "nd a better "tness value, but it will take more com-putational time. Also, there should be an optimal number of moves to obtain the optimal results for di!erent problems.

(3) &Cooling rate'

The &cooling rate' is the coe$cient to decide the rate of simulation &temperature' decrease. A slow &cooling rate' (e.g. 0)9) allows the molecules to align themselves into a completely ordered crystalline structure; this con"guration is the state of minimum energy for the system. If the &cooling' is too rapid (e.g. 0)1), the system does not reach the higher ordered state, but ends up in a high-energy state. The &cool-ing' schedule can be mathematically described as follows:

¹LCU"a¹MJB (15)

where: ¹LCU and ¹MJB are the simulation &temperatures' at the end and beginning of the &cooling' schedule; anda is the &cooling' rate, which can range from 0 to 1.

Fig. 2. Flowchart of simulated annealing

4.3. Representative design &chromosome'

The length of a design &chromosome' consists of a "xed number of binary digits. Also, the position and random number values in#uence the decoded value of the design

&chromosome'. To design a &chromosome' length to rep-resent an irrigation project, the cumulative numbers of crops within each command area are "rst calculated. Each crop is then assigned seven binary digits to repres-ent its area, which can range from 1 to 100%, in all of the

Fig. 3. A sample chromosome coding scheme to represent seven crops in the Delta project

percentage points, of the cumulative area in each com-mand area (seven binary digits give a value of 0}2!1, or 0}127 in decimal). Finally, the length of a design &chromosome' equals the cumulative number of crop types multiplied by seven.

For example, two command areas have been con-sidered in the Delta, Utah irrigation project for testing the model. The "rst command area, UCA1, includes three crop types and the second command area, UCA2, includes four crop types. Therefore, seven crop types are within these two command areas of the Delta, Utah irrigation project, and the length of a design some' should be 49. While considering a design &chromo-some' string of 49 binary digits, the seven crops in the two command areas can be depicted in coded form as shown in Fig. 3.

4.4. Decoding a design &chromosome' into a real number The design &chromosome' can be decoded into a deci-mal number to represent the crop area within each com-mand area. The conventional decoding method is used in this study. Consider a problem with k decision variables xG, i"1, 2,2, k, de"ned on the intervals xG 3[aG, bG]. Each decision variable can be decoded as a binary sub-string of length mG. The decoded decimal xG can be obtained from the following equation (McKinny & Lin, 1994): xG"aG#2bG!aGKG!1 KG QbQ;2Q (16) where: s is the summation identi"er for substring length. The following case study from the Delta, Utah project contains seven crop types in the two command areas. Therefore, this problem has seven decision variables xG, and i can range from 1 to 7. Without considering inherent crop area constraints, the percentage area of each crop type can range from 1 to 100% of the total command area. Therefore, the interval for each decision variable can be represented as xG 3 [1, 100], and aG equals 1 and bG equals 100. In conclusion, Eqn (16) can decode the binary digits into an actual number in the range from 1 to

100. The next step is to transfer this decimal number into crop area percentage AH%, and area AHF? within each

command area. A simple averaging technique was used, as given by

AH%" xH

,A HxH

100 (17)

AHF?"AH100%ASA? (18) where: j is the crop index; NA denotes the number of crops within each command area; and ASA? represents the area of each unit command area.

4.5. Rearrangement of design &chromosomes'

Rearrangement of design &chromosomes' is necessary to order the huge number of atoms within each simula-tion &temperature' value for minimizing the energy of the system. In this study, the rearrangements can be treated as a change in the location of the binary digits within the design &chromosome'; therefore, the newly allocated crop area can be obtained after decoding the rearranged de-sign &chromosome'.

There are three steps in the rearrangement of the design &chromosome': (1) move the binary digits at regions I}III depending on the length from the second cut site to the end of the design &chromosome'; (2) back the binary digits one position in region II; and (3) move the binary digits at regions III}I, depending on the length from the "rst digit to the "rst cut site. A design &chromosome' with 15 binary digits and two random break points demonstrates the procedure. The old and new design &chromosome', before and after rearrangement are shown in Fig. 4.

4.6. Annealing scheduling by boltzmann probability Annealing scheduling is the heart of the simulated annealing method. This procedure is the major di!erence from the traditional optimization methods (e.g. iterative improvement or Monte Carlo methods) that allows perturbations to move uphill in a controlled fashion;

Fig. 4. Rearrangement chromosome for the simulated annealing method

Fig. 5. Computer yow for annealing scheduling in the simulated annealing module

therefore, the simulated annealing method has the oppor-tunity to escape from a local optimum toward a global optimum. Figure 5 shows the computer #ow for anneal-ing schedulanneal-ing in the SA module and can be described in the following steps.

(1) The simulated annealing method allows many moves within one simulation &temperature' value; therefore, the "rst step is to compare the energy di!erenceDE (i.e. di!erence of project bene"t) from the previous move to the current move:

DE"EKMTC!EKMTC> (19) where: EKMTC and EKMTC> represent the project bene"t at current and previous moves.

(2) If the energy di!erence is negative (DE (0), the irrigation project maximum bene"t and related crop

areas are accepted at this move because the energy has been improved from the previous move to the current move.

(3) If the energy di!erence is positive (DE'0), this means the energy was not improved, but the irriga-tion project maximum bene"t and related crop areas still have the opportunity to update if the Boltzmann probability PP is greater than the generated uniform random number r. The Boltzmann probability can be de"ned as

PP"e\D#2Q? (20) From the above equation, it can be seen that PP is in#uenced by simulation &temperature' ¹Q?; that is, higher simulation &temperatures' will have higher PP values, and the system has a greater opportunity to update the con"guration ifDE'0. This also implies

Fig. 6. Seiver River Basin, Utah (Tzou, 1989). , Basin; , Rivers; , Cities

that the system at a higher simulation &temperature' has a higher ability to rearrange the atoms (i.e. to jump away from local optima) for "nding better, more optimal results. As the simulation &temperature' continues to decrease, the system tends to equilib-rium because the PP value is small, and there is no more ability to update the con"guration if*E '0. Finally, the global (or near global) optimum can be determined from this procedure.

4.7. Objective function and ,tness value

In this study, the objective function includes the in-come from crop harvest, cost of irrigation water and crop production cost. The objective is to maximize the irriga-tion project bene"t or the "tness value from the seven crops growing in the two command areas. Within the calculation loop of design &chromosome' size, the objec-tive function returns a "tness value to the model and then updates the "tness value and related crop-allocated area if this value is higher than previous ones. At the end of the design &chromosome' loop, the subsequent "tness value is the highest bene"t within the loop. Also, the maximum "tness value is selected from the generation number loop. Therefore, the "tness value and related crop area are the optimum results at the end of the calculations. The objec-tive function is mathematically expressed as

Maximize: , G ,A H(PGH>GH!SGH!FGH!OGH)AGH!= , G ,A HQGH (21) where: i, j is the command area and crop index; N is the number of command areas within irrigated project; NA is the number of crops within each command area; PGH is unit price of the jth crop in the ith command area in $ ha\; >GH is yields per hectare of the jth crop in the ith command area in ton ha\; SGH is seed cost per hectare of the jth crop in the ith command area in $ ha\; FGH is fertilizer cost of the jth crop in the ith command area in $ ha\; ¸GH is labour cost of the jth crop in the ith command area in $ ha\; OGH is operation cost of the jth crop in the ith command area in $ ha\; AGH is planted area of the jth crop in the ith command area in ha; = is unit price of irrigation water in $ m\; and QGH is cumu-lative water requirement of the jth crop in the ith com-mand area in m.

The objective function is subject to the following constraints.

(1) To consider social factors and to prevent one high -value crop from dominating the search for maximum

bene"t, the maximum and minimum area percentages must be considered for the crops:

A )A)A  (22)

where A  and A  are the minimum and maximum percentage area values of crop j in command area i in %, respectively.

(2) The cumulative water demand of crop j in command area i should be less than the available water supply for each command area:

,A

HQBCK)QQSN

(23) where: QBCK denotes the irrigation water requirement for crop j in command area i in m ; and QQSN repres-ents the available water supply for command area i in m.

Table 1

Recorded monthly weather data for the Delta, Utah

Month Temperature,3C No. of rainy days, day Rainfall,mm

Average SD Average SD Average SD

Jan !3)5 3)4 6)0 3)1 14)0 10)4 Feb 0)0 3)1 5)0 2)7 13)8 13)1 Mar 4)1 2)0 6)8 3)7 21)4 15)4 Apr 9)0 1)7 6)5 3)2 20)8 14)1 May 14)6 1)6 5)7 3)2 23)7 19)9 Jun 19)4 1)5 3)5 2)9 12)1 12)2 Jul 24)4 1)0 3)5 2)0 10)3 10)7 Aug 23)1 1)0 4)1 2)4 13)3 14)6 Sep 17)5 1)7 4)1 3)1 16)6 18)9 Oct 10)7 1)6 4)9 3)4 21)3 21)8 Nov 3)0 1)9 4)8 2)6 15)1 12)2 Dec !2)2 2)3 5)4 3)0 17)0 18)1

Fig. 7. The relationship between generated and recorded refer-ence crop evapotranspiration for Delta irrigated project in 1993:

, recorded; , generated

5. Application and results 5.1. Site description

The Wilson Canal System, close to the city of Delta in central Utah, was used in this study and is part of the many diversions in the Sevier River Basin operated by the Abraham Irrigation Company as an on-demand irri-gation system with a good communications network (Tzou, 1989). The Wilson Canal is 11 480 m in length with water being supplied from the Gunnison Bend Res-ervoir. Figure 6 depicts the location of the Sevier River Basin.

The climate in the Delta area is essentially a cold desert type, which is arid with cold winters and warm summers. The UCA1 and UCA2 command areas were selected within the Wilson Canal System for evaluating the model. The UCA1 command area has a 2896 m water course, 83)3 ha planted area, and three crop types are

planted: alfalfa, barley and maize. On the other hand, the UCA2 command area has a 12 350 m water course and 311)3 ha planted area. In addition, four crop types are planted: alfalfa, barley, maize and wheat.

5.2. Application of weather generation

Based on the monthly recorded weather data, the weather generation module can generate daily weather data for use in the on-farm irrigation scheduling module. The long-term recorded monthly meteorological data in Delta, Utah, as shown in Table 1, were used to generate the daily reference crop evapotranspiration, precipitation and mean, maximum and minimum air temperatures for use in the on-farm irrigation scheduling submodel. According to the test, the most suitable arid probability was 78 for the weather generation module generating the dried year and the seed number was 50 in this study. Figure 7 shows the relationships between the generated and 1993 recorded values of E2-. As presented in Fig. 7, the generated E2- was a little higher than the recorded E2- and has a similar trend.

5.3. Application to the on-farm irrigation scheduling module

The on-farm irrigation scheduling module deals with the daily simulation of the water balance to estimate irrigation water requirement and relative crop yield. Figure 8 shows the relationship between the daily soil moisture content, depth of irrigation and rainfall for alfalfa and barley crops in the UCA1 command area, respectively. Figure 9 compares the potential and cal-culated evapotranspiration for maize and wheat in the

Fig. 8. The relationship between soil moisture, irrigation depth and rainfall for (a) alfalfa and (b) barley crops in the UCA2 command area: , soil moisture; , irrigation depth

rainfall

Fig. 9. The relationship between potential ( ) and simulated ( ) evapotranspiration for (a) maize and (b) wheat in the

UCA4 command area

Table 2

Seasonal outputs for the UCA1 command area from the on-farm irrigation scheduling module; E2, evapotranspiration

Alfalfa Barley Maize

Potential E2, mm 1038)03 555)58 514)84

Actual E2, mm 907 505)71 460)94

Evaporation from wet soil surface, mm 2)08 21)35 13)37

Number of irrigations 6 4 3

Total irrigation depth, mm 1067)92 441)71 471)82

Deep percolation, mm 70)14 29)39 37)19

Surface runo!, mm 28)49 11)94 15)11

Yield reduction due to water stress, % 11)43 3)59 14)54

Yield reduction due to waterlogging, % 2)6 1)09 1)09

Relative crop yield, % 86)27 95)36 84)52

UCA2 command area. Tables 2 and 3 show the seasonal outputs from the on-farm irrigation scheduling module for the UCA1 and UCA2 command areas, respectively. The irrigation water requirements from the on-farm irri-gation scheduling module show that: (1) crops within the UCA1 command area are 1067)9, 441)7 and 471)8 mm for alfalfa, barley and maize, respectively; and (2) crops within UCA2 command area are 1039)5, 531)4, 490)9

and 539)4 mm for alfalfa, barley, maize and wheat, respectively. Also, the relative crop yield due to water stress and waterlogging can be summarized as follows: (1) crops within the UCA1 command area are 86)27, 95)36 and 84)52% for alfalfa, barley and maize, respectively; and (2) crops within UCA2 command area are 85)59, 95, 84)68 and 93)11% for alfalfa, barley, maize and wheat, respectively.

Table 3

Seasonal outputs for the UCA2 command area from the on-farm irrigation scheduling module; E2, evapotranspiration

Alfalfa Barley Maize Wheat

Potential E2, mm 1039)33 572)08 523)35 611)18

Actual E2, mm 906)15 528)68 469)61 558)05

Evaporation from wet soil surface, mm 3)38 37)85 21)88 34)4

Number of irrigations 7 6 4 6

Total irrigation depth, mm 1039)5 531)4 490)9 539)37

Deep percolation, mm 68)32 35)13 38)42 35)64

Surface runo!, mm 27)75 14)27 15)61 14)48

Yield reduction due to water stress, % 11)73 3)49 14)16 5)39

Yield reduction due to waterlogging, % 3)04 1)56 1)35 1)58

Relative crop yield, % 85)59 95)0 84)68 93)11

Table 4

Four data sets to test the simulated annealing submodel

Sets 1 2 3 4

Initial simulation &temperature' 1000 1000 1000 1000

Number of moves 50 70 90 110

&Cooling rate' 0)95 0)95 0)95 0)95

Number of runs 10 10 10 10

Table 5

Summarized results from four data sets for the simulated annealing method

Set 1 2 3 4 Standard deviation, $ 2103 1928 1581 2079 Average bene"t, $ 110 435 111 580 111 494 110 767 Maximum bene"t, $ 113 036 115 333 114 857 114 058 Minimum bene"t, $ 106 599 108 405 109 517 106 910 Number of moves 50 70 90 110

Initial simulation &temperature' 1000 1000 1000 1000

&Cooling rate' 0)95 0)95 0)95 0)95

5.4. Application of the simulated annealing

Three parameters are necessary for the simulated annealing method: (1) the initial simulation &temper-ature'; (2) the number of moves; and (3) the &cooling rate'. Marryott et al. (1993) stated that the total number of simulations required by simulated annealing is control-led by the length of the Markov chain (i.e. the number of moves). The Markov chain length represents the number of simulations per annealing temperature step required to ensure equilibrium in the optimization process. Therefore, Arts and van Laarhoven (1985) suggested that the minimum value of the Markov chain length (i.e. number of moves) must be satis"ed to ensure the location of an optimal (or near optimal) solution.

In this study, four rules were followed to "nd suitable parameters.

(1) The length of the Markov chain can be in the range of 10}100 times the number of decision variables. Also, the chain length can be cut or increased by a factor of 10 and ten simulated annealing runs performed for the given problem.

(2) The &cooling rate' can be chosen between 0)8 and 0)99 throughout the entire annealing run (Kirkpatrick et al., 1983; Kirkpatrick, 1984).

(3) A series of several runs can be performed with the same parameters because the result for each run is essentially di!erent and independent of the random starting point. If the standard deviation from all

Table 6

Simulated annealing results with initial simulation 9temperature: of 1000, number of moves 90, and 9cooling rate: of 0 '95 for the whole project and for two unit command areas (UCA)

Run Project UCA1 UCA2

Net benext, $

Water demand, ;103

m3

Alfalfa % Barley, % Maize, %

Water demand, ;103

m3

Alfalfa % Barley, % Maize, % Wheat %

Water demand, ;103 m3 1 112,715 3017)371 67)19 17)97 14)84 722)143 42)65 27)21 27)21 2)94 2295)228 2 114 857 3048)221 72)90 21)50 5)61 749)610 42)67 28)89 24)89 3)56 2298)612 3 109 517 2990)306 70)69 5)17 24)14 742)742 39)81 26)85 29)63 3)70 2247)564 4 110 593 2999)315 64)66 27)82 7)52 707)130 42)54 19)74 29)39 8)33 2292)185 5 112 255 3014)310 72)55 13)73 13)73 749)830 39)27 43)98 9)95 6)81 2264)481 6 110 005 3024)033 71)00 12)00 17)00 742)571 42)62 10)13 38)82 8)44 2281)462 7 112 646 3017)384 68)87 19)81 11)32 730)025 41)05 37)89 14)21 6)84 2287)359 8 110 144 2957)283 58)39 35)77 5)84 674)015 41)92 26)26 27)78 4)04 2283)267 9 110 318 2980)368 68)57 7)62 23)81 731)610 38)15 48)55 8)09 5)20 2248)758 10 111 886 2974)774 60)51 29)94 9)55 685)982 41)48 38)07 17)61 2)84 2288)793 Avg 111 494 3002'337 67'53 19'13 13'34 723'566 41'22 30'76 22'76 5'27 2278'771 Max 114 857 Min 109 517 SD 1581 S .-F . KUO E ¹ A ¸ .

Fig. 10. Sample graphs from the simulated annealing method to obtain optimal net benext with three parameters; initial temper-ature, TQ?"1000; cooling rate a"0)95; number of moves,_M"90 the runs is high, the parameters may not be suitable to "nd the optimal or near-optimal solutions. On the other hand, the parameters with low standard deviation from all the runs will have more con"dence to "nd the optimal or near-optimal results for the applied problems.

(4) Several sets of parameters are required to apply the model and the best set will have a higher average and a lower standard deviation of the project bene"t from a series of runs.

Based on the above rules, four data sets are used in this study to test the simulated annealing submodel as shown in Table 4. After ten runs were performed for each data set, Table 5 summarizes the "nal results from the four data sets. Table 5 indicates the average bene"t ranges from US$110 435 to 114 494 and the standard deviation ranges from US$2103 to 1581. For choosing the most suitable parameters in this study, set 3 is the best because it has the lowest standard deviation even if the average bene"t is a little lower than that for set 2. Therefore, the most appropriate parameters for the SA method in this study are as follows: (1) the initial simulation &temper-ature' equals 1000; (2) the number of moves equals 90; and (3) the &cooling rate' equals 0.95. Table 6 summarizes the "nal results from ten runs for these parameters. As presented, the maximum bene"t is up to US$ 114 857 on run 2, and standard deviation of bene"t is as low as US$ 1581 from the ten runs. Therefore, the parameters are sure to obtain the near global optimal values for this irrigation project planning problem.

Figure 10 displays the sample graph from the SA method to represent the bene"t from the Delta project during the searching process. Note that the graphs are read from right to left because the annealing proceeds from high to low simulation &temperatures'. At higher

simulation &temperatures', the annealing scheduling makes the graph climb either downward or upward, and the graphs tend toward equilibrium and reach the near global optimum results at the lower simulation &temperatures'.

6. Conclusion

This work develops a model based on the on-farm irrigation scheduling and the simulated annealing optim-ization method to provide guidelines on agricultural water resource planning and management. The model consists mainly of six basic modules: (1) a main module to direct the running of the model, with pull-down menu facility; (2) a data module to enter the required data by a user-friendly interface; (3) a weather generation module to generate the daily weather data; (4) an on-farm irriga-tion scheduling module to simulate the daily water requirement and relative crop yield; (5) a simulated annealing module to optimize the project maximum bene"ts; and (6) a results module to present results by tables, graphs and printouts. The model is applied to Delta, Utah for optimizing the maximum crop produc-tion bene"ts and identifying the crop area-allocated per-centages with the application of the simulated annealing method.

The simulated annealing method (SA) was imple-mented to optimize the agricultural water resource planning. As for the SA method, three parameters are needed in this method: (1) initial simulation &temper-ature'; (2) number of moves; and, (3) &cooling rate'. The SA module has been implemented with the on-farm irri-gation scheduling module to maximize the project bene-"t. The computational procedure of the SA can be divided into the following steps: (1) to receive the output from the irrigation module; (2) to enter the simulated annealing parameters through a user interface; (3) to de"ne the design &chromosome' representing the prob-lem; (4) to generate the random initial design &chromo-some'; (5) to decode the design &chromosome' into a real number; (6) to apply constraints; (7) to apply an objective function and a "tness value; (8) to implement the anneal-ing schedule by the Boltzmann probability; and (9) to set the &cooling rate' and criterion for termination.

The ability of annealing scheduling allows the SA method to update its con"gurations even though the energy is not improved from that of the previous simula-tion &temperature'. Therefore, the SA method can over-come the problem of traditional optimization methods that often get stuck in a local optimal because the traditional optimization method just allows the con"g-urations to be updated at the time of energy improve-ment from the previous iteration. In this study, the most

suitable data set has an average bene"t of US$ 114 500 and the standard deviation of bene"t equals US$ 1581. The relative controlling parameters are: (1) an initial simu-lation &temperature' of 1000; (2) the number of moves equal to 90; and, (3) a &cooling rate' of 0)95. The "nal results from the SA method for the Delta, Utah are: (1) a project bene"t of US$ 114 500; (2) project water demand of 3)002;10 m; (3) crop area percentages within the UCA1 command area of 67)5, 19)1, and 13)4% for alfalfa, barley and maize, respectively; and (4) crop area percentages within the UCA2 command area of 41)2, 30)8, 22)8, and 5)2% for alfalfa, barley, maize and wheat, respectively. References

Aarts E H L; Laarhoven P J M V (1985). Statistical cooling: a general approach to combinational optimisation problems. Philips Journal Research, 40, 193}226

Bohachevsky I O; Johnson M E; Stein M L (1986). Generated simulated annealing for function optimisation. Technomet-rics, 28(3), 209}217

Davis L (1991). Handbook of Genetic Algorithms. Van Nos-trand Reinhold, New York

Dougherty D E; Marryott R A (1991). Optimal groundwater management: 1. Simulated annealing. Water Resource Re-search, 27(10), 2493}2508

Hargreaves G L; Hargreaves G H; Riley J P (1985). Agricul-tural bene"ts for the Senegal River basin. Journal of Irriga-tion Drainage Engineering, ASCE, 111(3), 265}275

Hill R W; Keller A A; Boman B (1982). Crop yield models adapted to irrigation scheduling programs. Utah Agricul-tural Experiment Station, USU, Research Report 100

Keller A A (1987). The USU unit command area model. Water Management Synthesis II Project, WMS Report No. 71. USAID, Logan, UT

Kirkpatrick S (1984). Optimisation by simulated annealing: quantitative studies. Journal of Statistical Physics, 34(5), 975}986

Kirkpatrick S; Gelatt C D; Vecchi M P (1983). Optimisation by simulated annealing. Science, 220, 671}680

Laarhoven P J M V; Aarts E H L (1987). Simulated Anneal-ing: Theory and Applications. D. Rediel, Dordrecht, Holland

Maidment D R; Hutchinson P D (1983). Modeling water de-mands of irrigation project. Journal of Irrigation Drainage Engineering, ASCE, 109(4), 405}419

Marryott R A; Dougherty D E; Stollar R L (1993). Optimal groundwater management: 2. Application of simulated an-nealing to a "eld-scale contamination site. Water Resource Research, 29(4), 847}860

McKinny D C; Lin M D (1994). Genetic algorithm solution of groundwater management models. Water Resource Re-search, 30(6), 1897}1906

Otten R H J M; Ginneken L P P P V (1989). The Annealing Algorithm. Kluwer Academic, Norwell, MA

Prajamwong S (1994). Command area decision support system for irrigation projects. PhD Thesis, Utah State University, Logan, UT

Smith M (1991). CROPWAT: Manual and Guidelines. FAO of UN, Rome, Italy

Tzou C T (1989). Veri"cation and calibration of the USU hydraulic model. MS Dissertation, Utah State University, Logan, UT

Walker J C (1992). Simulated annealing applied to the peanut growth model for optimisation of irrigation

schedul-ing. PhD Thesis, North Carolina State University,