Stochastic Programming in Water Resources System Planning: A Case Study and a Comparison of Solution Techniques

Working Paper

S T O C W C PROGRAMMING IN WATER RESOURCES

s!BTEM

PLANNING:

A

CASE

STUDY

AND A

COMPARISON OF SOLUTION TECHNIQUES

J. DupaEoviL A. G a i v o r o n s k i 2. Kos

T.

S z b n t a i

August 1986 WP-86-40

International Institute for Applied Systems Analysis

A-2361 Laxenburg, Austria

NOT FOR QUOTATION WITHOUT THE PERMISSION OF THE AUTHORS

STOCHASTIC PROGR4MMING IN WATER RESOURCES

SYSTEM

PLANNING:

A

CASE

STUDY

AND

A

COMPARISON OF SOLUTION TECHNIQUES

J. DupaEoviL A. Gaivoronski 2. Kos

T. Szbntai

August 1986 WP-06-40

Working Papers are interim r e p o r t s on work of t h e International Institute f o r Applied Systems Analysis and have received only limited review. Views or opinions e x p r e s s e d h e r e i n do not necessarily r e p r e s e n t those of t h e Institute o r of i t s National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS 2361 Laxenburg, Austria

FOREWORD

Stochastic programming methodology is applied in this p a p e r to a c a p i t a l in- vestment problem in water r e s o u r c e s . The a u t h o r s introduce possible model formu- lations and t h e n find t h e solutions by using a number of specific solution tech- niques. These are p a r t l y based on t h e SDS/ADO t a p e f o r s t o c h a s t i c programming problems and also on s t a n d a r d l i n e a r and nonlinear programming packages. Such a n a p p r o a c h allows a thorough analysis of t h e solution as well as a comparison of t h e algorithmic p r o c e d u r e s used to obtain t h e s e solutions.

Alexander B. Kurzhanski Chairman System and Decision Sciences Program

AUTHORS

Dr. J i t k a ~ u p a C o v 6 , r e s e a r c h s c h o l a r at SDS in 1985-1986 i s p r o f e s s o r at t h e Department of S t a t i s t i c s , Faculty of Mathematics and Physics, Charles University, P r a g u e .

Dr. Alexei Gaivoronski, r e s e a r c h s c h o l a r at SDS i s s e n i o r r e s e a r c h e r of t h e V.

Glushnov Institute of Cybernetics, Kiev, USSR.

Ing. Dr. ~ d e n z k K o s , r e s e a r c h s c h o l a r at IIASA i n 1980, i s s e n i o r r e s e a r c h e r a t t h e Faculty of Engineering of t h e Technical University, Prague.

Dr. Tomds Szdntai i s p r o f e s s o r a t t h e Department of Operations R e s e a r c h of t h e Eotvos Lordnd University, Budapest.

ABSTRACT

To analyze t h e influence of increasing needs upon a given water r e s o u r c e s system in E a s t e r n Slovakia and to g e t a decision on t h e system development and ex- tension, s e v e r a l s t o c h a s t i c programming m o d e l s c a n b e used. The t w o s e l e c t e d models are based on individual probabilistic c o n s t r a i n t s f o r t h e minimum s t o r a g e and f o r t h e f r e e b o a r d volume supplemented by one joint probabilistic c o n s t r a i n t on releases or by a nonseparable penalty t e r m in t h e objective function. Suitable numerical techniques f o r t h e i r solution are applied to a l t e r n a t i v e design parame- ter values. As a r e s u l t , t h e p a p e r gives a n answer to t h e case study which i s based on multi modeling within t h e framework of s t o c h a s t i c programming and, at t h e s a m e time, i t gives a comparison of various solution techniques p a r t l y included in t h e SDS/ADO collection of s t o c h a s t i c programming codes.

-

vii

-

STOCHASTIC PROGRAMMING IN WATER RESOURCES

SWlXM

PLANNING:

n

CASE

SIVDY AND A

COWPARISON OF SOLUTION TECHNIQUES

J. h p a E ovk A. Gaivoronski, 2. Kos a n d T. S z h n t a i

1. mRODUCTION

Operation of r e s e r v o i r s in water r e s o u r c e s system i s a multistage s t o c h a s t i c control problem. Water r e s o u r c e s planning h a s to consider multiple u s e r s and ob- jectives, r e s e r v o i r operation policies need t o b e analyzed to obtain a n effective use of water r e s o u r c e s . In t h e planning p r o c e d u r e , a g r e a t v a r i e t y of water r e s o u r c e s systems designs or operation plans need t o be confronted and t h e i r economic, environmental and social impacts evaluated. Two basic t a s k s are often solved: (a) a new system i s developed o r an existing one i s enlarged by s o m e pro- posed investments in r e s e r v o i r s , pipes, canals, pumping stations, hydroelectric water plants, etc., or (b) t h e management and c o n t r o l of a n existing system h a s t o be a l t e r e d t o accommodate t o t h e new conditions. Analysis in both t h e s e c a s e s rests on mathematical modeling, t h e objectives and constraint of t h e problem have t o b e e x p r e s s e d mathematically. The mathematical models involve t h e selection of many engineering, design and operating variables. The optimization means t h e determi- nation of t h e best values of t h e s e variables regarding t h e constraints. A lot of t h e variables and p a r a m e t e r s o c c u r r i n g in t h e objective function and in t h e con- s t r a i n t s has t o b e taken as s t o c h a s t i c as they d e s c r i b e a s t o c h a s t i c real life prob- l e m .

Out of many possible goals, t h e water supply f o r industry and irrigation, flood control and r e c r e a t i o n purposes a r e considered in this study. The m o s t important decision variable is t h e s t o r a g e capacity of r e s e r v o i r s .

T h e r e were drawn up a lot of stochastic programming models f o r inventory control and water s t o r a g e problems by P r d k o p a (1973). One of t h e s e models w a s applied f o r designing s e r i a l l y linked r e s e r v o i r systems (Prdkopa e t a1 1978).

Another application of t h e s e models concerning t h e flood c o n t r o l r e s e r v o i r design

w a s described by P r d k o p a and Szdntai (1978). In Czechoslovakia where t h e analyzed system is located, t h e r e i s a n old tradition t o use c h a n c e constrained models f o r t h e considered t y p e of problems. Applications of t h e s e models were described by DupaEov6 and Kos (1979), Kos (T979) in conjunction with t h e l i n e a r decision r u l e o r with t h e d i r e c t control. Here, we shall follow t h i s tradition and in addition, we shall consider t h e possible use of a l t e r n a t i v e stochastic programming models.

I t i s t h e evident s t o c h a s t i c n a t u r e of inflows which causes t h e most of t h e sta- tistical and modeling problems and which causes t h e necessity of a detailed analysis of t h e d a t a and of t h e variables occurring in any of model formulations.

In water r e s o u r c e s modeling, f o u r types of variables o c c u r :

(1) constant coefficients and p a r a m e t e r s t h a t are used f o r design values (e.g., active s t o r a g e of t h e r e s e r v o i r , reliability, c o s t coefficients) o r variables with small variation (e.g., t h e withdrawal of water f o r t h e thermal power sta- tion);

(2) uncontrolled random v a r i a b l e s with a known o r estimated distribution (month- ly flows, meteorologic variables, e.g., precipitation, potential evaporation);

(3) p a r t l y controlable random variables with known distribution (e.g. irrigation water requirements with controlable a c r e a g e ) ;

(4) random variables with a n incomplete knowledge of distribution, such as t h e fu- t u r e demands, f u t u r e p r i c e s and costs.

In o u r p a p e r t h e type (3) of stochastic variables were incorporated into t h e model taking into account t h e i r joint probability distribution. In this way t h e in- t e r r e l a t i o n s between t h e succeeding months values of irrigation water require- ments could be considered. I t is t h e t h e o r y of logconcave measures developed by Prdkopa (1971) which gives t h e t h e o r e t i c a l background f o r t h e handling of joint probability c o n s t r a i n t s in s t o c h a s t i c programming problems.

W e pursued in this study two objectives. One i s t o develop optimization sto- chastic models t o determine t h e r e q u i r e d r e s e r v o i r capacity under increasing needs. The goal of t h e model i s t o show how t h e irrigation, flood c o n t r o l and r e - c r e a t i o n needs influence t h e r e s e r v o i r operation and t h e reliability in meeting t h e multiple goals of t h e water r e s o u r c e s system. This was done in section 2 , 3 and 6 which contain t h e description of the water r e s o u r c e s system, a s well a s s e v e r a l suitable s t o c h a s t i c programming models and a n application of one of them t o t h e basin of the Bodrog River in t h e E a s t e r n p a r t of Slovakia (Czechoslovakia), s e e

Figures 1 and 2.

Another of o u r objective was t o use this comparatively simple but meaningful and r e a l life problem t o t e s t various a p p r o a c h e s f o r solution of s t o c h a s t i c pro- gramming problems and coh,.pare various solution techniques which a r e implement- ed in IIASA and constitute SDS/ADO collection of s t o c h a s t i c programming codes.

The r e s u l t s of t h i s comparison c a n be found in sections 4 and 5.

2. THE WATER RESOURCES SYSTEM

The water r e s o u r c e s system consists of t h r e e r e s e r v o i r s , V., D., K., two of which are in operation (D., V.) and t h e t h i r d one is to be built o r not (K.)

-

s e e Fig- u r e 3. The main purposes of t h e water r e s o u r c e s system are t h e irrigation water supply, industrial uses

-

^mainly water withdrawal f o r t h e thermal power station, flood control ( b e t t e r flood alleviation), environmental conservation and r e c r e a - tion. The subsystem of t h e Laborec r i v e r (see Figure 4) w a s originally designed f o r industrial water supply and flood control. During t h e operation of t h e r e s e r v o i r V., a n a c c e l e r a t e d development of t h e r e c r e a t i o n o c c u r r e d and t h e demands f o r maintenance of t h e minimum r e c r e a t i o n pool during t h e summer period were sup- ported by authorities. The area of irrigation grew and according to t h e plan will be increased substantially. The main questions f o r t h e decision-makers are:

-

Can the p r e s e n t e d water r e s o u r c e s system still m e e t a l l t h e requirements?

And if s o , with which reliability?

-

Is t h e construction of t h e r e s e r v o i r K. n e c e s s a r y and when it will b e neces- s a r y ?

The analysis of this problem w a s divided into two s t e p s . The f i r s t one comprises t h e screening modeling and i t is discussed in this p a p e r . Optimization models, however, cannot r e f l e c t all t h e details of t h e water r e s o u r c e s system operation. T h e r e f o r e t h e r e s u l t s of t h e stochastic optimization model a r e supposed t o be verified using t h e stochastic simulation model with t h e input g e n e r a t e d by t h e methods of stochastic hydrology. For t h e stochastic programming screening model, which i s t h e s u b j e c t of t h e p r e s e n t p a p e r , a n aggregated model w a s used and t h e monthly flows and t h e i r r i g a t i o n water requirements were aggregated into f o u r periods:

Prague

C Z E C H O S L O V A K I A 0

v-

Brno

FIGURE

1 Location of the water resources system of the Bodrog river.

POLAND

HUNGARY

FIGURE 2 The water r e s o u r c e system allocation.

(1) November till April of t h e following y e a r , (2) May and June,

(3) July and August,

(4) September and October.

The f i r s t period starts at t h e beginning of t h e hydrological y e a r and comprises t h e winter and t h e s p r i n g periods filling t h e r e s e r v o i r s . As t h e main i n t e r e s t of t h e i r - rigation i s c o n c e n t r a t e d on t h e vegetation period and r e c r e a t i o n season, t h e aggregation of t h e six months is acceptable. The second and f o u r t h periods in- clude i r r i g a t i o n and industrial demands, t h e t h i r d period includes in addition t h e r e c r e a t i o n demands. The requirements for t h e minimum pool due to environmental

Industry

Thermal Power Plant

v

Flood Control

FIGURE 3 Schematic r e p r e s e n t a t i o n of t h e water r e s o u r c e s system of t h e Bodrog r i v e r .

c o n t r o l a n d enhancement a n d flood c o n t r o l p e r t a i n to a l l t h e p e r i o d s .

3. THE MATHEMATICAL

YODELS

The models were designed f o r s c r e e n i n g a l t e r n a t i v e s on t h e cost basis.

Using probabilistic c o n s t r a i n t s , t h e r e s u l t identifies t h e c a p a c i t y z o of r e s e r v o i r V. t h a t should meet t h e needs with a p r e s c r i b e d reliability. The t a s k i s formulated as t h e cost of r e s e r v o i r V. minimization. As t h e cost i s a n i n c r e a s i n g

/

A

\ Reservoir K.

/

\

Laborec R.

Reservoir V . (Recreation)

Flood Control

FIGURE 4 Schematic r e p r e s e n t a t i o n of t h e subsystem of t h e Laborec r i v e r .

function of r e s e r v o i r c a p a c i t y zo, w e c a n evidently minimize t h e c a p a c i t y zo in- s t e a d of minimizing t h e cost.

The c o n s t r a i n t s involve exceeding of water supply need

pi +

^di by t h e r e l e a s e d volume zi in periods 2, 3 and 4 (vegetation periods). The needs consist of t h e fixed demand d i (minimum flow and industrial water needs) and random demand

pi

(irrigation water requirements). A s to t h e f i r s t period, t h e fixed demand d l

(caused mainly by t h e needs of t h e thermal power station) should be met with such a high probability t h a t t h e deterministic constraint zl 2 d l w a s used. Taking into account t h e i n t e r c o r r e l a t i o n s of random demands

p i ,

ⁱ

=

2 , 3, 4, t h e constraint f o r t h e vegetation period, w a s formulated as follows

where a is t h e r e q u i r e d joint probability level.

The second type of c o n s t r a i n t s r e f l e c t s t h e environmental c o n t r o l and enhancement, fishing and r e c r e a t i o n needs. These requirements are e x p r e s s e d in t h e form of maintaining a minimum pool o r minimum r e s e r v o i r s t o r a g e m i . In periods 1 , 2 and 4 , t h e environmental and fishing needs a r e r e f l e c t e d , in period 3 t h e r e c r e a t i o n needs are added. This constrain1 i s t h e n as follows

where sf i s t h e r e s e r v o i r s t o r a g e , ai t h e r e q u i r e d individual probability t h r e s - holds in period i

.

The t h i r d t y p e of c o n s t r a i n t e x p r e s s e s t h e flood control t a r g e t assuming t h a t some flood control s t o r a g e v i will b e f r e e in t h e r e s e r v o i r operation during t h e whole period i with probability y i , i.e.,

The resulting optimization problem h a s t h e form minimize zo

and s u b j e c t t o additional c o n s t r a i n t s

which stem mostly from n a t u r a l hydrological and morphological situation. The u p p e r bounds u f , i

=

1,

. . .

, 4 , are given by t h e volume of flood with r e s p e c t i v e duration and probability of exceedance.

Using t h e d i r e c t (zero-order) decision r u l e and neglecting losses due to eva- poration, t h e r e s e r v o i r s t o r a g e s sf c a n b e e x p r e s s e d via t h e water inflows and r e l e a s e s in t h e r e l e v a n t periods. Let r j denote t h e water inflow in t h e j - t h period and l e t

ti

denote t h e cumulated water inflow,

Denote f u r t h e r by s o t h e initial r e s e r v o i r s t o r a g e at t h e beginning of t h e hydro- logical y e a r . A s a r u l e , w e c a n p u t s o

=

m 4, i.e., t h e r e s e r v o i r s t o r a g e i s supposed t o b e at i t s minimum a f t e r t h e vegetation period. Repeated use of t h e continuity equation gives

Substituting into (2) and (3) yields

Using t h e corresponding lOOpX quantiles z k ( p ) of t h e distribution of t h e random v a r i a b l e s t k , k

=

1 , 2, 3, 4 , w e c a n r e w r i t e t h e individual probabilistic c o n s t r a i n t s (7) and (8) in t h e form

Unfortunately, t h i s simple d e v i c e d o e s not a p p l y to t h e joint p r o b a b i l i s t i c con- s t r a i n t (1).

The resulting optimization problem minimize z o

c a n b e solved e.g. by s p e c i a l t e c h n i q u e s developed by P r d k o p a et al. (1978), see S e c t i o n 4.

Alternatively, s t o c h a s t i c programming decision models c a n b e built solely o n t h e evaluation a n d minimization of t h e o v e r a l l e x p e c t e d c o s t s , which contain not only t h e cost of t h e r e s e r v o i r of t h e c a p a c i t y s o b u t a l s o t h e losses c o n n e c t e d with t h e f a c t t h a t t h e n e e d s were not fulfilled a n d / o r t h e r e q u i r e m e n t s on t h e minimum r e s e r v o i r s t o r a g e a n d o n t h e flood c o n t r o l s t o r a g e were not met. This t y p e of models i s called mostly stochastic programs w i t h recourse or w i t h penalties.

Suppose f i r s t t h a t t h e c o n s t r a i n t s o n w a t e r supply, minimum w a t e r s t o r a g e a n d flood c o n t r o l s t o r a g e d o not c o n t a i n random v a r i a b l e s , i.e., w e h a v e (besides of t h e u p p e r a n d lower bounds (5))

In case of random p i , i

=

2 , 3, 4 a n d

t k ,

k

=

1,

. . .

^, 4, t h e c h o s e n decision z o , z

. . .

, z need n o t fulfill c o n s t r a i n t s (12) f o r t h e a c t u a l ( o b s e r v e d ) values of p i ,

t k .

If t h i s i s t h e c a s e , costs evaluating losses o n c r o p s (not being i r r i g a t e d o n a sufficiently high level), o n t h e d e c r e a s e of r e c r e a t i o n (due to t h e lower r e s e r v o i r pool) a n d t h e economic losses d u e to flood c a n b e a t t a c h e d to t h e d i s c r e p a n c i e s .

L e t c ( z o ) d e n o t e t h e cost of r e s e r v o i r of t h e c a p a c i t y z o , let t h e penalty functions b e of t h e t y p e

2 0 a n d n o n d e c r e a s i n g if y

>

⁰

.

Denote by v:, v 2 ,

vz

t h e p e n a l t y functions c o r r e s p o n d i n g to t h e c o n s i d e r e d t h r e e t y p e s of c o n s t r a i n t s in (12), i

=

2, 3 , 4 , k

=

1,

. . .

^, 4. W e try to find s u c h a deci- sion f o r which t h e t o t a l e x p e c t e d cost will b e minimal s u b j e c t to inequality con- s t r a i n t s (5):

4

minimize c ( z o )

+

E _{v,'(@i + di}

-

zi )

s u b j e c t t o i o S z o 5 U O

The c h o i c e of t h e p e n a l t y functions should b e b a s e d on a d e e p economic a n d environmental a n a l y s i s of t h e underlying problem. On t h e o t h e r hand, f o r a s c r e e n i n g s t u d y , i t s e e m s s a t i s f a c t o r y to r e s t r i c t t h e c h o i c e to piece-wise l i n e a r or piece-wise q u a d r a t i c p e n a l t y functions.

a ) piece-wise l i n e a r p e n a l t y (simple r e c o u r s e model) All penalty functions 9 are of t h e form

~ ( y )

=

qy', where q 2 0 a n d y'

=

max (0, y )

.

The coefficient q h a s t o b e given by t h e decision maker; i t c o r r e s p o n d s t o t h e unique costs f o r t h e d i f f e r e n t c o n s i d e r e d d i s c r e p a n c i e s . A s a r e s u l t , w e h a v e t o

s u b j e c t t o t h e c o n s t r a i n t s (5).

The solution method i s mostly based on approximation of t h e marginal d i s t r i - butions of

P i , t

by d i s c r e t e o n e s a n d , i n case of c ( z o ) l i n e a r , t h e resulting p r o - g r a m c a n b e solved by simplex method with u p p e r bounded v a r i a b l e s ( s e e e.g. t h e SPORT p r o g r a m (by N a z a r e t h ) contained in t h e ADO/SDS t a p e ) . F o r o t h e r t y p e s of piece-wise l i n e a r p e n a l t i e s see e.g. DupaEov& (1980).

b ) If t h e c o s t function c ( z o ) i s s t r i c t l y q u a d r a t i c on t h e c o n s i d e r e d i n t e r v a l

<lo, u O > , a s p e c i a l t y p e of piece-wise q u a d r a t i c p e n a l t y function was suggested by R o c k a f e l l a r a n d Wets (1985), namely,

~ ( y ) = O f o r y S O

= - y 2 / p 1

2 ^{f o r 0 5 y 5 p q}

( s e e Figure 5).

F o r solving problem (15) with p e n a l t y functions of t h e mentioned t y p e , p r o - g r a m LFGM (by King) implementing R o c k a f e l l a r a n d Wets' Lagrangian Finite Gen- e r a t i o n Method i s at disposal on t h e ADO/SDS t a p e . The p a r a m e t e r s p , q of t h e in- dividual penalty functions h a v e t o b e given by t h e decision m a k e r .

The use of t h e s t o c h a s t i c p r o g r a m (15) with p e n a l t i e s d o e s n o t r e q u i r e t h e knowledge of t h e joint d i s t r i b u t i o n of t h e random n e e d s f o r i r r i g a t i o n water. I t means, t h a t d u e t o t h e assumed s e p a r a b i l i t y of t h e penalty functions (i.e., d u e t o t h e f a c t t h a t s h o r t a g e in i r r i g a t i o n w a t e r i s penalized in e a c h of t h e t h r e e vegeta- tion p e r i o d s s e p a r a t e l y a n d t h e t o t a l penalty is t a k e n as t h e sum o v e r t h e t h r e e

FIGURE 5

p e r i o d s ) , n o i n t e r c o r r e l a t i o n s are considered. Alternatively, we c a n a t t a c h a penalty cost to t h e s i t u a t i o n , when t h e t o t a l needs f o r water i r r i g a t i o n in t h e vege- tation p e r i o d as a whole w e r e not met. In t h a t c a s e , w e c a n t a k e e.g.

with p of t h e form (12) a n d minimize

s u b j e c t to c o n s t r a i n t s (5).

Finally, i t i s possible to combine t h e p r o b a l i s t i c c o n s t r a i n t s a n d t h e penaliza- tion: o n e c a n define t h e set of admissible decisions by means of p r o b a b i l i s t i c con- s t r a i n t s (4) and inequalities (5) a n d , at t h e same time penalize t h e o c c u r r e n c e of t h e d i s c r e p a n c i e s in (12) by c o r r e s p o n d i n g penalty terms in t h e o b j e c t i v e function.

In t h i s case s t u d y , t h e set of admissible solutions w a s defined by t h e individual p r o b a b i l i s t i c c o n s t r a i n t s on a minimum pool a n d o n t h e flood c o n t r o l s t o r a g e , i.e., by t h e system of inequalities (9), ( l o ) , a n d by inequalities (5). Instead of t h e joint p r o b a b i l i s t i c c o n s t r a i n t on t h e water supply needs, o n e penalty t e r m of t h e form

i ( B i +

d i

-

x 2 . i

=

2 . 3 . 4 )

=

c [ m a r ( B i + d i - x i ) + ]

i = 2 , 3 , 4

0, max

( B f +

d i

-

^{X i )}

i = 2 , 3 , 4

w a s used. The resulting problem

minimize x o

+

_i = 2 , 3 , 4 max

( B i +

d i

-

^Z*

1 +}

k

s u b j e c t t o

z

^{xi S z k ( l}

^-

^{a k )}

⁺

^{m 4}

^-

^{m k , k}

⁼

^{1 ,}

^{. . .}

^{, 4 ,}

i = l

c a n be solved by t h e s t o c h a s t i c quasigradient method Ermoliev ( 1 9 7 6 ) o r by tech- niques designed t o solution of t h e complete r e c o u r s e problem; s e e Section 4 . For t h e optimal solution of ( 1 7 ) , t h e values of t h e joint probability

were computed.

Observe t h a t in ( 1 7 ) , t h e c o s t c ( x o ) of t h e r e s e r v o i r of capacity z o is sup- posed t o be l i n e a r in t h e i n t e r v a l L o S z 0 S u o and t h a t t h e coefficient c of t h e penalty t e r m evaluates unit losses due to water s h o r t a g e r e l a t i v e t o t h e c o s t p e r unit capacity of t h e r e s e r v o i r .

4. SOLUTION TECHNIQUES USED TO SOLYE THE

PROBLEM

In this section we shall d e s c r i b e briefly the solution techniques used f o r nu- merical experiments with models ( 1 1 ) and ( 1 7 ) described in t h e previous section.

We could choose among s e v e r a l stochastic optimization programs from t h e SDS/ADO stochastic optimization l i b r a r y available at IIASA.

F o r solving problem (11) with t h e joint p r o b a b i l i t y c o n s t r a i n t , nonlinear p r o - gramming t e c h n i q u e s c a n b e used. The c h o i c e among them d e p e n d s on t h e p r o p e r - t i e s of t h e set of f e a s i b l e solutions. F o r log-concave p r o b a b i l i t y m e a s u r e s , which i s t h e case of multidimensional normal, gamma, un:form, D i r i c h l e t d i s t r i b u t i o n s of

8,

t h e set d e s c r i b e d b y

i s convex and among o t h e r s , t h e method of f e a s i b l e d i r e c t i o n s , s u p p o r t i n g h y p e r - plane method a n d p e n a l t y methods supplemented b y a n e f f i c i e n t r o u t i n e f o r comput- ing t h e values of t h e function p ( z ) a n d i t s d e r i v a t i v e s c a n b e applied. (For a s u r - vey see P r d k o p a (1978), (1986).) F o r multinormal d i s t r i b u t i o n , t h e s u p p o r t i n g hy- p e r p l a n e method was implemented at IIASA by Szdntai as PCSP c o d e .

F o r solving problem ( 1 7 ) , o n e possibility i s to a p p r o x i m a t e t h e o r i g i n a l d i s t r i - bution of

8

^by a s e q u e n c e of discrete distributions. The penalty c a n b e given im- plicitly via so c a l l e d second stage program: f o r f i x e d values of z , ,

pi,

d , , i

=

2 , 3 , 4 ,

? ( p i +

d ,

-

x i , i

=

2 , 3 , 4 ) e q u a l s to t h e optimal value of t h e o b j e c t i v e func- tion in t h e following l i n e a r p r o g r a m

minimize c y

so t h a t t h e problem (17) c a n b e c o n s i d e r e d as t h e complete r e c o u r s e problem a n d solved accordingly. F o r d i s c r e t e d i s t r i b u t i o n t h e a p p r o a c h l e a d s to l i n e a r pro- gramming p r o b l e m s of s p e c i a l s t r u c t u r e which should b e e x p l o i t e d by a d e q u a t e solution techniques. One s u c h t e c h n i q u e i s L-shaped algorithm by Van S l y k e a n d Wets (1969) implemented at IIASA by Birge in NDSP code. F o r f u r t h e r exposition see e.g. Kall (1979), Wets (1983).

Finally, to s t o c h a s t i c programming models of e x p e c t a t i o n t y p e s u c h as (17).

s t o c h a s t i c q u a s i g r a d i e n t (SQG) method c a n b e applied, see Ermoliev (1976), Ermo- l i e v a n d Gaivoronski (1984). The s t o c h a s t i c optimization s o l v e r ST0 based on t h i s method s o l v e s t h e following problem

minimize E, f ( z , w )

=

F ( z )

s u b j e c t to c o n s t r a i n t s

where f ( z , o ) i s a functic? which depends on decision v a r i a b l e s z a n d random p a r a m e t e r s o. The u s e r h a s to p r o v i d e a n algorithmic d e s c r i p t i o n of t h i s function.

The set X i s t h e set of c o n s t r a i n t s a n d in t h e c u r r e n t implementation, i t would b e t h e set defined by l i n e a r c o n s t r a i n t s , s a y ,

The optimal solution of (18)-(19) i s r e a c h e d i t e r a t i v e l y s t a r t i n g from a n initial point z 0 by applying t h e following i t e r a t i v e p r o c e d u r e :

w h e r e p, i s t h e s t e p s i z e ,

CS -

s t e p d i r e c t i o n ,

9 -

p r o j e c t i o n o p e r a t o r o n t h e set X:

The p r o j e c t i o n in t h e ST0 i s p e r f o r m e d using QPSOL q u a d r a t i c programming pack- a g e , see Gill et a l . (1983).

The s t e p d i r e c t i o n

CS

should, roughly s p e a k i n g , b e in a v e r a g e c l o s e to t h e g r a d i e n t of t h e o b j e c t i v e function F ( z )

=

E,j'(z, o ) at point z S , although indivi- dual

CS

may b e f a r from a c t u a l values of t h e g r a d i e n t . This i s e x p r e s s e d with t h e h e l p of conditional e x p e c t a t i o n s :

where a, i s some vanishing t e r m . Each p a r t i c u l a r s t r a t e g y of choosing s e q u e n c e of s t e p s i z e s p, a n d s t e p d i r e c t i o n s l e a d to p a r t i c u l a r algorithm a n d many s u c h s t r a t e g i e s are implemented in t h e p r o g r a m ST0 some of them are f a i r l y sophisticat- e d . I t i s a l s o possible to c h a n g e s t r a t e g i e s i n t e r a c t i v e l y during optimization pro- c e s s . Detailed d e s c r i p t i o n of t h i s p r o g r a m i s given in Ermoliev a n d Gaivoronski (1984), f o r t h e o r e t i c a l b a c k g r o u n d a n d f u r t h e r r e f e r e n c e s see Ermoliev (1976).

H e r e w e s h a l l d e s c r i b e only f e a t u r e s r e l e v a n t to t h e numerical e x p e r i m e n t s con- d u c t e d with water r e s o u r c e models.

F i r s t of all i t w a s n e c e s s a r y t o put t h e problem in t h e f o r m (18)-(19). Observe t h a t t h e model with t h e joint probabilistic c o n s t r a i n t cannot b e easily put in t h e frame-work (18)-(19) while t h e expectation models (14), (15). (16) and (17) are al- r e a d y fc-mulated in t h e r e q u i r e d fashion. But expectation models d o not give t h e value of probability

which i s important f o r decision making. T h e r e f o r e w e conducted numerical e x p e r i - ment in t w o s t a g e s .

1 Solve t h e problem based on t h e expectation model:

minimize

E f l ( z , 8 ) = z o + c E m a x max

[pi

+ d i - z i j

i = 2 , 3 , 4

I

s u b j e c t to c o n s t r a i n t s

and additional c o n s t r a i n t s (5) (see (17)). Program ST0 g e n e r a t e s c e r t a i n ap- proximation z to

*

t h e optimal solution.

2 Evaluate t h e value of probability c o n s t r a i n t , t h a t i s compute

This w a s done using t h e g e n e r a t o r of multivariate normal distribution from IMSL l i b r a r y which w a s a l s o used on t h e f i r s t s t a g e . If p ( z * ) a p p e a r e d to b e less than admissible level a then w e i n c r e a s e d t h e penalty coefficient c from (21) and solved t h e problem (21)-(22) again. This p r o c e s s w a s r e p e a t e d until d e s i r a b l e value of ( z

*

) w a s r e a c h e d . Clearly, p ( z * ) ⁴ 1 when c w.

N o w w e s h a l l discuss s t r a t e g i e s which w e r e used f o r solving (21)-(22). In t h i s case i t i s possible t o compute subgradients j',(zS, p S ) of function f ( z , 8 ) f o r given values of z

=

z S and random p a r a m e t e r s

8 = pS,

t h e components are

and t a k e

tS

= f , ( z S ,

PS)

in method ( 2 0 ) . W e decided, however, to use a method which i s based only on values of f ( z S ,

PS)

and i s applicable t h e r e f o r e f o r m o r e complex models. Various t y p e s of random s e a r c h techniques and finite d i f f e r e n c e s are implemented in STO. For t h i s model t h e following analog of random s e a r c h w a s used:

-

G e n e r a t e L random v e c t o r s ?fl,

. . .

^, tlSL : tlSi

=

(adi,.

.

^,

a

f ' ) where tl;i are uniformly and independently distributed on [- 6 , , 6 , ] .

-

G e n e r a t e 1 independent random v e c t o r s p l ,

. . .

,

pSL, psi =

( @ i f ,

@ti)

with multivariate normal distribution.

-

Compute

p :

f z i ( z s ,

B S ) =

During computations w e took 1 ^S L ^S 5 . - c if @;

+

d i - 2 ;

=

=

maxto, max

18; +

d i - z ; j j i

0 otherwise

-

S e l e c t p, and p e r f o r m o n e s t e p in ( 2 0 ) .

The value of s e a r c h s t e p w a s t a k e n constant 6 ,

=

1 0 . Stepsize p, w a s updated interactively, but a f t e r s e v e r a l t r i a l s t h e following p a t t e r n emerged, which w a s re- peated a f t e r w a r d s :

1 S s S 5 0 5 0 S s S 1 5 0 1 a f t e r w a r d s

In o r d e r to g e t e x a c t solution i t is n e c e s s a r y to t a k e 6 , ^--, 0 . But in t h i s particu- l a r case i t a p p e a r e d t h a t sufficiently good approximation w a s obtained even with fixed 6 , . Usually i t took 200-300 i t e r a t i o n s to g e t solution. The r e s u l t s obtained by t h i s method for solving water r e s o u r c e problems with d i f f e r e n t p a r a m e t e r s are discussed f r o m t h e application point of view in t h e section 6. One of t h e sets of problem p a r a m e t e r s w a s used to compare performance of various models and solu- tion techniques d e s c r i b e d at t h e beginning of this section.

5. COMPARISON OF THE SOLUTION TECHNIQUES

This i s t h e p a r t i c u l a r problem of t h e type (17) which w a s used f o r comparison:

minimize F ( z )

= ^EJ

( z , o )

=

zo

+

^{C E ,}

[ [

max 0 , max ^1=2,3,4 toi

+

^di

-zii

s u b j e c t to c o n s t r a i n t s

where di

=

12.7, i

=

2, 3, 4. The random v e c t o r o

=

( o l , o z , u3) i s distributed nor- mally with

expectations e

=

(20.2, 27.37, 10.65) s t a n d a r d deviations q

=

(8.61. 10.65. 6.00)

0.360 0.125 and c o r r e l a t i o n matrix 0.360 1. 0.571 ,

[::I25 0.571 1.

1

penalty coefficient c was equal 100. This problem i s t h e same as discussed in t h e case study in Section 6 , only t h e u p p e r bound u o

=

500 mil. m3 was used instead of u o

=

334 mil. m 3 and t h e reliabilities T k

=

0.75 were t a k e n instead of y k

=

0.4.

The convenient f e a t u r e of t h i s problem i s t h a t we c a n easily obtain v e r y good lower bound f o r solution by minimizing zo s u b j e c t to c o n s t r a i n t s s t a t e d above. This gives F ( z * ) 2 494.88 where z * i s t h e optimal point.

5.1. A p p r o z i m a t i o n s of s t o c h a s t i c problem b y Large s c a l e Linear p r o g r a m m i n g problem

To use this a p p r o a c h i t i s necessary to approximate initial continuous distri- bution by d i s c r e t e distribution consisting of N points. Various ways of doing this are described in Birge, Wets (1986). H e r e w e used t h e following two schemes:

1 "Intelligently"

-

variance matrix w a s computed from t h e given c o r r e l a t i o n matrix and variances;

-

eigenvalues r i and eigenvectors v i , i

=

1, 2 , 3 of t h e variance matrix were computed;

-

number k w a s chosen and f o r one-dimensional normal distribution points bi , i

=

1,

. . .

^{, k} were s e l e c t e d such t h a t

f o r i

=

^1,

...

^k

-

¹

-

approximating d i s c r e t e distribution consisted of all points

f o r a l l i = l , .

. .

, k , j = l , .

. .

, k , 1 = l , .

. .

, k and

These points were t a k e n with equal probability l / ( k 3 ) . Approximating distri- bution chosen in t h i s way i s symmetric and h a s t h e same expectation and vari- a n c e matrix as original distribution.

2 J u s t throwing specific number N of points distributed according t o original distribution and assign t o e a c h of them probability 1 / N .

After such discretization i s performed, t h e original problem becomes equivalent t o t h e following stochastic optimization problem with complete r e c o u r s e :

C N minimize z o

= - ^vj

j = l

s u b j e c t to c o n s t r a i n t s

and c o n s t r a i n t s (24).

where uj E j

=

1,

. . .

,

N ,

are t h e points where t h e d i s c r e t e distribution i s concentrated.

This i s l i n e a r programming problem which could b e of v e r y high dimension be- cause l a r g e N i s needed f o r a c c u r a t e approximation. W e did not, however, make v e r y accurate approximations and t h e r e f o r e were a b l e to use a g e n e r a l purpose l i n e a r programming tool, namely l i n e a r programming p a r t of MINOS 4.0, see Mur- tagh and S a u n d e r s (1983).

The r e s u l t s are summarized in t h e Table 1. To e a c h e n t r y in t h e l e f t column of this t a b l e c o r r e s p o n d t w o rows of numbers in t h e columns 2-4. The u p p e r number c o r r e s p o n d s to "intelligent" approximation and t h e lower to "just throwing". The column marked ST0 p r e s e n t s r e s u l t s obtained by t h e v a r i a n t of SQG d e s c r i b e d in section 4 a f t e r 100 i t e r a t i o n s , e a c h i t e r a t i o n r e q u i r e d 40 observations of random function. I t i s amazing t h a t f o r much bigger problem of 1331 points in case of "in- telligent" approximation MINOS found solution much f a s t e r t h a n f o r 343 points. I t i s a l s o worth noting what a big difference makes "intelligent" approximation as com- p a r e d with "just throwing": approximation which u s e s only 125 points is b e t t e r than one with 1331 points. The estimate of t h e value of t h e objective function in Table 1 w a s made using a sample of 10000 points g e n e r a t e d by subroutine f o r multivariate normal distribution from t h e IMSL l i b r a r y .

W e complete t h e discussion by comparison of t h e empirical expectation and c o r r e l a t i o n matrix computed using t h e s a m e 10000 random numbers with t h e corresponding preasigned p a r a m e t e r values. I t gives a n idea how a c c u r a t e t h e es- timates a r e :

expectations of random v e c t o r o e

=

^{(20.2000, 27.3700, 10.6500)}

empirical expectations from 10000 points (20.2855, 27.3912, 10.6841) 0.360 0.125

'

c o r r e l a t i o n matrix of o 0.360 1. 0.571 1 2 0.571 1.

,

empirical c o r r e l a t i o n matrix

1. 0.354429 0.107631 0.354429 1. 0.572717 .0.107631 0.572717 1.

1

TABLE 1

- -

Total n u m b e r of points

in approximating d i s t r i b u t i o n

Optimal value of t h e complete r e c o u r s e o b j e c t i v e function Estimate of t h e value of o r i g i n a l o b j e c t i v e function Number of MINOS i t e r a t i o n s MINOS u s e r time

5.2. Stochastic q u a s i g r a d i e n t method

In t h e ST0 implementation t h e u s e r c a n c h o o s e between two options: i n t e r a c - t i v e a n d automatic. In t h e i n t e r a c t i v e option t h e u s e r c h o o s e s t h e s t e p s i z e p, from ( 2 0 ) himself, h i s judgment i s b a s e d on s u c h notions as "oscillatory b e h a v i o r " of v a r i a b l e s o r "visible t r e n d " a i d e d by some additional "measures of p r o c e s s s quali- ty". These are displayed on t h e terminal along with c u r r e n t point and u s e r c a n in- t e r r u p t i t e r a t i v e p r o c e s s s a n d c h a n g e t h e value of s t e p s i z e . In o r d e r to use t h i s possibility e f f e c t i v e l y t h e u s e r h a s to b e q u i t e e x p e r i e n c e d . However, as p r a c t i c e shows e v e n a n i n e x p e r i e n c e d u s e r c a n quickly g e t n e c e s s a r y skills. All t h i s i s d e s c r i b e d in more d e t a i l in Gaivoronski (1986).

In e x p e r i m e n t s with i n t e r a c t i v e mode t h e s t e p d i r e c t i o n from ( 2 0 ) was com- p u t e d using c e r t a i n a n a l o g u e of random s e a r c h : g e n e r a t e random v e c t o r h

=

( h l , h 2 , h 3 ) w h e r e hi a r e independent random v a r i a b l e s uniformly d i s t r i b u t e d on [ - G , G ] a n d G i s c h o s e n i n t e r a c t i v e l y from t h e i n t e r v a l [I, 201, compute llh

I(

^corn-

p u t e v e c t o r u

=

( u l , up, u g ) s u c h t h a t

where oi a r e independent o b s e r v a t i o n s of o

r e p e a t t h e f i r s t two s t e p s L times with d i f f e r e n t independent o b s e r v a t i o n s of r a n - dom v a r i a b l e s w a n d h , o b t a i n 1 v e c t o r s u and t a k e as s t e p d i r e c t i o n s

t s

a v e r a g e of t h e s e v e c t o r s (in example given below 1

=

10).

S t e p s i z e w a s c h o s e n i n t e r a c t i v e l y . A t f i r s t sufficiently l a r g e s t e p s i z e c h o s e n (10.0) which w a s r e d u c e d if i r r e g u l a r b e h a v i o r of t h e z v a r i a b l e s w a s o b s e r v e d . The r e s u l t s a r e in Table 2.

TABLE 2

S t e p S t e p s i z e zo z 1 2 2 2 3 z 4 F(z)

number

The value of s t e p G in random s e a r c h w a s s e t 10.0

A t t h i s point t h e s t e p in random s e a r c h w a s c h a n g e d t o 5.0

During a c t u a l computations information was displayed more f r e q u e n t l y and, of c o u r s e , without l a s t column which c o n t a i n s t h e e s t i m a t e s of t h e values of t h e objec- tive function f (z) at t h e c u r r e n t points using t h e same 10000 o b s e r v a t i o n s of t h e

random v e c t o r o which were used f o r estimations of MINOS r e s u l t s . These estimates were obtained a f t e r w a r d s (and consumed a lot of computer time!). From t h e table i t i s evident t h a t v e r y good solution ( b e t t e r than 1331 points approximate scheme) w a s obtained a f t e r 20 i t e r a t i o n s , t h a t is a f t e r 800 function evaluations and a l l sub- sequent points oscillated in close vicinity.

The major disadvantage of t h e i n t e r a c t i v e option is t h a t i t r e q u i r e s too much from t h e u s e r . T h e r e f o r e automatic option w a s developed in which computer simu- l a t e s behavior of a n e x p e r i e n c e d u s e r . F o r experiment with automatic option w e choose t h e simplest version of SQG: o n e which uses f o r s t e p d i r e c t i o n t h e values of

I,

( z , o ) , which are v e r y e a s y to compute h e r e : f * ( z , o )

=

I l . , o., - c t i , - c t 2 , - c t 3 j where t j

=

1 if

O j + d j t i - z j + I =maxIO, _{i =2,3,4} max I o i - l + d i - z i j j and t j

=

0 otherwise, c

=

100

-

penalty coefficient.

On e a c h i t e r a t i o n only o n e observation of

I,

( z S , o ) w a s computed, which w a s used as s t e p direction.

Stepsize ps w a s computed automatically according to t h e simple rule:

-

on e a c h i t e r a t i o n o n e observation

I

( z S , oS ) was made and t h e s e observations were used to compute estimate F ( s ) of t h e c u r r e n t value of t h e objective function:

c u r r e n t p a t h length L ( s ) w a s computed also:

S

L ( s )

= C

^{J J z i}

+'

^-zil(

i =1

-

initial stepsize p i w a s chosen sufficiently l a r g e (in examples below i t w a s 5.0) and e a c h M i t e r a t i o n s t h e condition f o r reducing stepsize w a s c h e c k e d (in ex- amples below M

=

20, t h a t is conditions were c h e c k e d on i t e r a t i o n s number 20, 40, 60,.

. .

). This condition i s t h e following:

p s + l = D p s if ( F ( s - K ) - F ( s ) ) / ( L ( s ) - L ( S -K)) S A PS + I

=

PS otherwise

In examples below D

=

0.5, K

=

^{20, A}

=

0.01, s t a r t i n g point was (1000, 100, 100, 100, 100). Each i t e r a t i o n r e q u i r e d one observation of random function f ( z , 0 ) and one observation of i t s gradient. Below t h e r e are two r u n s with dif-

f e r e n t sequences of random v e c t o r s o. One i s "very good" a n o t h e r

-

^"not as good but quite reasonable". When c u r r e n t point a p p r o a c h e s optimum t h e e v e n t It,

=

1 j becomes less and l e s s likely. T h e r e f o r e t h e method spends much of i t e r a t i o n s standing at t h e same point ( f o r instance i t e r a t i o n s 260-400 in Tables 3 and 4). The l a s t column again w a s obtained a f t e r w a r d s using t h e same 1 0 0 0 0 observations of w as in t h e previous tests. I t shows t h a t t h e algorithm r e a c h e s quite good vicinity of solution a f t e r 120 i t e r a t i o n s (except a slight jump on i t e r a t i o n 420). In Table 4 w e have a big jump on i t e r a t i o n number 220, t h e method, however, quickly r e a c h e s t h e vicinity of solution again. This jump is due to too big stepsize. After i t e r a t i o n number 240 everything is OK again.

The problem with SQG i s t h e stopping c r i t e r i o n and c u r r e n t l y t h e c r i t e r i o n implemented i s based o n assumption t h a t if stepsize becomes'too small w e are in t h e vicinity of optimum. This of c o u r s e i s not n e c e s s a r i l y t r u e but e x p e r i e n c e shows t h a t n e v e r t h e l e s s t h i s c r i t e r i o n is quite r e l i a b l e if coupled with r e p e a t e d runs.

The same input d a t a were used f o r solving t h e problem (11) with joint proba- bility constraint. N e w v a r i a b l e s S o

=

^{z o}

-

^{L o and}

=

^{z 1}

-

d l were introduced t o eliminate t h e individual lower bounds. For t h e resulting program

minimize S o

O S ~ o S ~ O - ~ O ; O S ~ l S ~ ~ - d ~ ; O S ~ ~ S ~ ~ ; O S ~ ~ S ~ ~ ; O S ~ ~ S ~ ~ ,

t w o solution techniques were implemented:

s u b j e c t to

P

z Z S 0 1

+

d 2 z 3 2 0 2 + d 3

2 o3

+

d 4

2 a

TABLE 3 "Very good".

S t e p Stepsize 20 Z 1 Z~ 2 3 z 4 F(z

1

number

TABLE 4 "Not as good as previous but quite reasonable".

S t e p Stepsize

Z o

^{Z 1} 2 2 3 4 F(z

1

number

5.3. The PCSPcode b y SzdLntai

The p r o g r a m s o l v e s problems of s t o c h a s t i c programming with joint p r o b a b i l i t y c o n s t r a i n t s u n d e r assumption of multinormal distribution of random right-hand s i d e s ( o f

= pi,

i

=

2, 3 , 4 in o u r c a s e ) . I t i s based on Veinott's s u p p o r t i n g i . 1 ~ - p e r p l a n e algorithm ( s e e Veinott (1967)). The individual u p p e r bounds on v a r i a b l e s a r e handled s e p a r a t e l y a n d t h e p a r a m e t e r s of t h e multinormal distribution are used t o g e t a s t a r t i n g f e a s i b l e i n t e r i o r solution. F o r c o n s t r u c t i n g t h e n e c e s s a r y l i n e a r a n d s t o c h a s t i c d a t a f i l e s , o n e c a n t u r n t o t h e brief documentation by Ed- wards (1985). Some computational r e s u l t s of t h e test problem with d a t a given at t h e beginning of t h i s S e c t i o n are given in Table 5.

TABLE 5 Results of t h e calculations by using t h e PCSP code.

zo 1

Zz

3 4 P r o b . lev. CPU time No. of

cutting planes

5.4. The a p p l i c a t i o n of t h e n o n l i n e a r v e r s i o n of t h e MINOS s y s t e m

F o r t h i s p u r p o s e o n e h a s t o write a s e p a r a t e s u b r o u t i n e named CALCON which c a l c u l a t e s t h e value of t h e p r o b a b i l i s t i c c o n s t r a i n t s a n d i t s g r a d i e n t . The s u b r o u - t i n e s are contained in t h e PCSP code. They w e r e coded on t h e b a s e of a n improved simulation technique by Szdntai (1985). Computational r e s u l t s given in t h i s s e c t i o n show t h a t t h e d i r e c t a p p l i c a t i o n of t h e MINOS system f o r t h e solution of t h e optim- ization problem (11) i s l e s s comfortable t h a n t h e use of t h e PCSP code. We obtained t h a t without giving a good initial s e t t i n g on t h e values of t h e nonlinear v a r i a b l e s t h e MINOS system failed t o find a f e a s i b l e solution.

Finally in Table 7 t h e r e i s t h e summary of experiments.

In t h e f i r s t column t h e r e i s a s h o r t d e s c r i p t i o n of e x p e r i m e n t including name of t h e p r o g r a m , number of approximating points a n d t y p e of approximation ( f o r MINOS, I a n d R mean t h e same as b e f o r e ) , number of i t e r a t i o n s a n d indication of in- t e r a c t i v e or automatic mode ( f o r STO), value of p r o b a b i l i t y c o n s t r a i n t ( f o r PCSP).

TABLE 6 Results of t h e calculations by using t h e nonlinear MINOS system d i r e c t - ly. (We applied t h e z 2

=

63.035, z3 =77.345, z 4

=

30.0 initial s e t t i n g s . F o r t h e s e initial values t h e p r o b a b i l i s t i c c o n s t r a i n t is infeasible with value 0.866.)

2 3 z 4 P r o b . lev. CPU time No. of major i t e r a t i o n s

In t h e column 7 t h e r e are values of o b j e c t i v e function (24) computed by a v e r a g i n g 1 0 000 o b s e r v a t i o n s g e n e r a t e d by IMSL s u b r o u t i n e f o r multivated normal distribu- tion. In t h e column 8 are t h e values of probability c o n s t r a i n t computed using t h e same random numbers a n d in t h e l a s t column c p u time of t h e VAX 780. F o r i n t e r a c - t i v e mode of ST0 c p u time i s not included b e c a u s e t h e c r u c i a l f a c t o r t h e r e w a s u s e r ' s r e s p o n s e .

Experiments s u g g e s t t h a t both c o n s i d e r e d models g i v e c o m p a r a b l e r e s u l t s and t h a t t h e methods ST0 a n d PCSP contained in t h e ADO/SDS l i b r a r y of s t o c h a s t i c p r o - gramming c o d e s p e r f o r m e d b e t t e r on t h i s p a r t i c u l a r problem t h a n t h e d i r e c t u s e of s t a n d a r d LP o r NLP p a c k a g e s . To s o l v e t h e c a s e s t u d y , t h e s t o c h a s t i c q u a s i g r a - d i e n t method w a s c h o s e n as i t s implementation i s n o t essentially limited by t h e as- sumed t y p e of d i s t r i b u t i o n a n d t h r o u g h i n c r e a s i n g t h e penalty coefficient value a n a p p r o x i m a t e optimal solution which fulfills t h e joint p r o b a b i l i t y c o n s t r a i n t c a n b e achieved.

6.

THE CASE STUDY

In t h e c a s e s t u d y of t h e water r e s o u r c e s system in t h e Bodrog R i v e r basin t h e model (17) w a s used with a d d e d c o n s t r a i n t s ( s e e discussion):

The input values of t h e 3-dimensional multinorrnal distribution of were as follows (with e x c e p t i o n of c o r r e l a t i o n matrix in mil. m3)

TABLE 7 Experiment

MINOS MINOS MINOS MINOS MINOS MINOS ST0 I ST0 I ST0 A 1 ST0 A2 PCSP PCSP PCSP PCSP PCSP PCSP PCSP

Function value

Prob.

value

CPU time

C o r r e l a t i o n m a t r i x

The p a r a m e t e r s of t h e m a r g i n a l normal d i s t r i b u t i o n s of c u m u l a t e d monthly inflows

< i w e r e

As t h e v a l u e s a k

=

O.g, k

=

1, 2, 3, 4 a n d y k

=

0.4, k

=

1, 2, 3, 4 w e r e used ( s e e d i s c u s s i o n ) t h e following v a l u e s of q u a n t i l e s zk w e r e o b t a i n e d (in mil. m3

T h e v a l u e s di w e r e as follows (mil. m3) P e r i o d k

zk (1.

-

^0.9)

zk (0+4>

The minimum a n d maximum r e s e r v o i r c a p a c i t y z o was: l o

=

100 mil. m3 a n d u,,

=

³³⁴ mil. m3. (In t h e a l t e r n a t i v e d i s c u s s e d in t h e p r e v i o u s s e c t i o n a less r e a l i s t i c v a l u e u O

=

⁵⁰⁰ mil. m3 was u s e d . ) The u p p e r bounds f o r v a r i a b l e s w e r e 252 mil. m3 ( u i = 252, i

=

1 , 2, 3, 4 in mil. m3) t h a t i s t h e volume of a long-term flood. This c o n s t r a i n t was n o t e f f e c t i v e a n d t h e r e f o r e i t was n o t a n a l y s e d .

1 2 3 4

146.8 205.0 252.9 283.0 272.5 342.1 397.1 446.1

Due to r e c r e a t i o n p u r p o s e s , t h e a c c e p t a b l e minimum s t o r a g e in t h e 3-rd p e r i o d i s m 3

=

1 9 4 mil. m3. However, t h e comparison of t h e t h i r d inequality of (9), t h e t h i r d inequality of (10) t o g e t h e r with z o S 334 g i v e s a n u p p e r b o u n d of 1 8 9 , 4 fc; t h e sum m 3

+

v3, so t h a t t h e p a r a m e t e r v a l u e m 3

=

1 9 4 would l e a d to c o n t r a d - i c t o r y c o n s t r a i n t s . T h a t s why t h e minimum s t o r a g e v a l u e m 3 h a s b e e n p u t up to 1 3 7 mil. m3 ( s e e a l t e r n a t i v e C) a n d t h e s t o r a g e v a l u e s m k , k

=

1 , 2, 3, 4 h a v e b e e n k e p t f i x e d o v e r a l l p e r i o d s ( s e e a l t e r n a t i v e s A a n d B). The r e l i a b i l i t y of maintaining t h e summer r e s e r v o i r pool h a s b e e n e v a l u a t e d e x p o s t .

6.1. C h o i c e of r e l i a b i l i t y v a l u e s

The v e r y i m p o r t a n t p a r a m e t e r s of t h e model are t h e r e q u i r e d p r o b a b i l i t i e s a , 7 2 a n d a i . The v a l u e a i s t h e r e q u i r e d joint p r o b a b i l i t y of t h e water s u p p l y . The tests with t h e model h a v e shown t h a t i t i s n e c e s s a r y to a d d t h e d e t e r m i n i s t i c con- s t r a i n t (25) in o r d e r to s e c u r e t h e r e q u i r e d v a l u e s of c o n s t a n t i n d u s t r i a l water demands. The p e n a l t y t e r m i s o f t e n so weak t h a t t h e c o n s t r a i n t (25) may b e violat- ed. Using t h e d e t e r m i n i s t i c c o n s t r a i n t (25) a r e l a t i v e l y low v a l u e a, e.g. a

=

0.85 may b e a c c e p t a b l e . However, t h e v a l u e of joint p r o b a b i l i t y (1) i s t h e o u t p u t of t h e model; t h e r e f o r e t h e condition (1) i s t e s t e d a n d if i t c a n n o t b e fulfilled, t h e a l t e r - n a t i v e i s r e j e c t e d .

The v a l u e s a i r e f e r to t h e r e l a t i v e l y s t r o n g environmental a n d t e c h n i c a l re- q u i r e m e n t s f o r maintaining t h e minimum r e s e r v o i r pool. T h e r e f o r e a f = 0.9, i

=

1, 2 , 3, 4 w a s c h o s e n .

The c h o i c e of t h e v a l u e s 7 * w a s r a t h e r difficult. They r e f e r to t h e i m p o r t a n t c o n s t r a i n t s imposed o n t h e r e s e r v o i r V o p e r a t i o n t h a t a r i s e from flood c o n t r o l re- q u i r e m e n t s t h a t s t i p u l a t e t h a t a c e r t a i n s p a c e

-

flood c o n t r o l s t o r a g e , b e held empty. This r e q u i r e m e n t c a n n o t b e e a s i l y e x p r e s s e d i n t h e model d u e to i t s a g g r e - g a t e d c h a r a c t e r . The p r o b a b i l i t y t h a t t h e f r e e b o a r d s t o r a g e i s empty means a l s o t h a t t h e r e i s n o s p i l l during t h i s p e r i o d . If t h e p e r i o d i s s h o r t , e.g. o n e d a y , t h e p r o b a b i l i t y 7 may c o r r e s p o n d to t h e r e q u i r e d r e l i a b i l i t y . F o r l o n g e r p e r i o d s e . g . o n e month, t w o months, half a y e a r , e x t r e m e l y l a r g e s t o r a g e c a p a c i t i e s will b e n e c e s s a r y if n o s p i l l s h a l l t a k e p l a c e d u r i n g t h i s p e r i o d with p r o b a b i l i t y e q u a l 0.75. (This v a l u e h a s b e e n used in t h e p r e v i o u s s e c t i o n a n d t h e t o t a l r e s e r v o i r c a p a c i t y z o as high as a p p r o x . 5 0 0 mil. m3 w a s n e c e s s a r y . ) T h e r e f o r e t h e flood c o n t r o l p r o b l e m s are o f t e n t r e a t e d in a s e p a r a t e model a n d t h e r e q u i r e d p r o b a b i l i -

t i e s 7 a r e a d a p t e d to t h e r e s u l t i n g v a l u e s of t h i s s e p a r a t e model. I t i s n e c e s s a r y t o i n t e r p r e t p r o p e r l y t h e meaning of t h e s e p r o b a b i l i t i e s . F o r i n s t a n c e t h e p r o b a b i l i t y 7

=

0.4 of maintaining t h e f r e e b o a r d s t o r a g e d o e s n o t s a y t h a t t h e flood c o n t r o l i s on a low level b u t t h a t t h e f r e e b o a r d s p a c e will b e filled (and possibly some spill may o c c u r ) with t h i s p r o b a b i l i t y during t h e p e r i o d c h o s e n . With t h i s f a c t in mind and a c c o r d i n g to t h e r e s u l t s of t h e s e p a r a t e flood c o n t r o l model t h e value yi

=

0.4 i

=

1, 2 , 3, 4 w a s c h o s e n in t h i s s e c t i o n .

6.2. R e s u l t s

In t h e f i r s t r e s u l t i n g a l t e r n a t i v e A of t h e r e s e r v o i r V design a n d o p e r a t i o n t h e i n p u t p a r a m e t e r s were:

mk = 5 7 mil. m3, k

=

1, 2 , 3, 4

-

minimum r e s e r v o i r s t o r a g e ,

v k

=

7 0 mil. m3, k

=

1 , 2 , 3, 4 - flood c o n t r o l s t o r a g e ( f r e e b o a r d s t o r a g e ) ,

which g a v e

x o

=

291.6 mil. m3

-

t h e t o t a l r e s e r v o i r c a p a c i t y x I

=

107.9 mil. m3

-

t h e t o t a l r e l e a s e in t h e 1-st p e r i o d x 2

=

6 9 . 6 mil. m3

-

t h e t o t a l r e l e a s e in t h e 2-nd p e r i o d x 3 = 69.8 mil. m3

-

t h e t o t a l r e l e a s e in t h e 3-rd p e r i o d x 4

=

35.7 mil. m3

-

t h e t o t a l r e l e a s e in t h e 4-th p e r i o d .

S t a r t e d at s o

=

m 4

=

5 7 mil. m3, using t h e computed optimal r e l e a s e s and t h e cumulated monthly inflows e q u a l to t h e q u a n t i l e s zk (1

-

0.9) and zk (0.4) we c a n compute t h e c o r r e s p o n d i n g s t o r a g e s si (in mil. m3):

P e r i o d i

I

^{1 2 3}

The underlined v a l u e s are c r i t i c a l , i . e . , t h e c o r r e s p o n d i n g inequalities in (9) and (10) a r e fulfilled as e q u a t i o n s . A s t h e existing t o t a l r e s e r v o i r c a p a c i t y i s 304 mil.

m3, we o b s e r v e from t h e a c t i v e c o n s t r a i n t s

I

min s t o r a g e s i m a x s t o r a g e s i

5 7 9 5 . 9 84.5 62.6 57.0 5 7 221.6 221.6 206.8 219.9