QUANTILE OPTIMIZATION TECHNIQUES WITH APPLICATION TO CHANCE CONSTRAINED PROBLEM FOR WATER-SUPPLY SYSTEM DESIGN Andrei V.Naumov Department of Industrial and Operations Engineering University of Michigan Ann Arbor, Mi 48109 Andrei I.Kibzun Department of Applied Mathematics Moscow Aviation Institute Moscow, 127080, USSR Technical Report 92-5 January, 1992. 1

QUANTILE OPTIMIZATION TECHNIQUES WITH APPLICATION TO CHANCE CONSTRAINED PROBLEM FOR WATER-SUPPLY SYSTEM DESIGN ANDREI I.KIBZUN AND ANDREI V.NAUMOV Department of Applied Mathematics, Moscow Aviation Institute, Moscow, Russia, 127080 Department of Industrial and Operations Engineering, University of Michigan, Ann Arbor, Michigan 48109 ( Revised on December 26, 1991.) This paper considers an optimimization problem in water supply systems, which can be formulated as a minimization with probabilistic constraints. A nontraditional approach for the solution of this problem using quantile optimization algorithms is presented. In this case the General Minimax Approach allows to use the linear program solution techniques for obtaining an upper bound on the optimal value of the objective function. Numerical results and application of the stocastic quasi gradient algorithm for asymptotic exact solution are also discussed. ( WATER-SUPPLY SYSTEM DESIGN; PROBABILISTIC CONSTRAINTS; GUARANTEEING SOLUTION; QUANTILE OPTIMIZATION ALGORITHMS) 1. Introduction. Optimization in water-supply systems often arises in the design of distant solar powered plants for fresh water production when local fresh water sources are absent. 1

The model considers stochastic deviations of the solar plant power and the possibility of monthly external fresh water supplies. The objective function of this problem is the sum of the capital cost of solar powered plant construction and the cost of external water supplies. We consider the probability of complete satisfaction of the demands to be restricted. This optimization problem is the typical stochastic programming problem with chance constrained. Other similar applications where considered for example in Prekopa (1980), Prekopa and Ganczer, Deak, Patyi (1980), Prekopa and Szantai, Prekopa and Kelle (1978). In general, methods for the solution of such problems are based ( see Prekopa 1987) on the application of suitable nonlinear programming techniques supplied by Monte Carlo simulation procedures to find values and gradients of the probability functions. We refer to a few papers ( Prekopa and Szantai, Prekopa 1987, Prekopa and Kelle 1978, Mayer 1980, Mayer 1988), where the interested reder can find descriptions of how to use nonlinear optimization techniques for such problems. In our paper we reduce the initial chance constrained problem to an equivalent unconditional quantile function optimization, i.e. optimization of an auxiliary random variable distribution quantile function without additional probabilistic restriction. This allows us to use the General Minimax Approach (GMA) (see Malyshev Kibzun 1987) to use high developing linear program solution technique to obtain the guaranteeing solution, i.e. the upper bound on the optimal value of the objective function. 2

This may be practically useful because the computational requirements for the Monte Carlo simulation procedures depends greatly on the probabilistic constraints confidence level a. In this problem, a is interpreted as the reliability of the system, and so a is close to unity. Computational requirements for the guaranteeing solution procedure do not depend on a, moreover the upper bound tends to an exact optimal value of the objective function as a increases (see Malyshev and Kibzun 1987). We also discuss the application of the stochastic quasi gradient algorithms for the quantile function minimization (see Kibzun 1991, Kibzun and Kurbakovskiy 1991a, Kibzun and Kurbakovskiy 1991b) to find asymptotically with probability one an exact solution of the equivalent problem. Numerical results are considered in the last section. 2. Problem statement. Let h = (h,..., hn) be a nonrandom vector of monthly water demands over n months and let u' = (u1,..., u,) be a nonrandom vector of monthly external supply plans. The internal water-supply system consists of a solar powered plant for water production with area S and cistern of volume V, where the remaining water is kept. The solar powered plant's productivity depends on monthly solar activity and can be expressed as a random vector w = (w1,..., wn). The probability distribution of random vector w is known. The total cost of the water-supply system is given by: n S = io(S, V, u') = aS + a2V ao E uj (1) j=1 where al is the cost per m2 of the solar powered plant; 3

a2 is the cost per m3 of the cistern; ao is the cost per m3 of externally delivered water. The surplus water at the end of jth month is expressed by: xj = min(xj_l, V) + Swj + u - h, xo = O, j 1, n, (2) where min(xj_-, V) is water remaining from the previous month, which can not exceed the volume of the cistern. The system works normally if the following condition is satisfied: Xj >0, j = Ln. (3) Since equation (2) includes the random vector w, then satisfying (3) is random event. The optimization problem is to find nonnegative system parameters S, V and nonnegative external water supply parameter u, which minimize the objective function (1). The optimal parameters are: u* = arg min 0o(i), (4) iEU by given values of ao, al, a2, h and distribution w, subject to a probabilistic constraint: P(u) = P{w: min xj(i, w) > 0} > a. (5) 1,n Here a E (0, 1) is the given system reliability and u EU = {(u1,..., un, S, V): S > O,V >O, > O, j l,n}. Using (2) and(3) we can transform (5) to the probabilistic constraint of simultaneously satisfying of N = n * (n + 1)/2 bilinear inequalities, because every jth month 4

we must add j inequalities to the previous inequality system. Let us explain this fact in an example over n=2 months. Then (5) is: x1 = wlS + ul - hi > 0 P w: >a. (6) x2 = min(xl, V) + w2S + 2 - h2 >_ 0 In fact, this system of inequalities is equivalent to the probability of satisfying N = 3 inequalities of the following form: w1S + u1 - hi > 0 P w: (w + w2)S + ul + u2- h - h2 > 0 a (7) V + W2S + 2 - h2 > 0 Hence problem (4),(5) is a stochastic program with a linear objective function and probabilistic constraints of satisfying N bilinear inequalities. 3. Reduction to a quantile optimization problem. Let us reduce the problem (4),(5) to an equivalent unconstrained quantile optimization problem, i.e optimization problem of a confidence level a quantile function for some auxiliary random variable distribution, without additional probabilistic constraints. Let us express variable ul using (S,V,o) and (uj, j = 2,n) from (1). Define also a new variable z = V - 4&/2a2. After changing variables, the probability constraint (5) can be reduced to the following form: P{w: oi(u, w) <, i= 1,N + 2} > a, (8) 5

where u E U = (S,z, u2,...,un):S > O,uj > O, j = 2,n}, and functions (i have the structure: <.i(.) = Sfi(w) + ]i(z, U2..., ), U = 1, N IN+1(') = -2z (9) (N+2(-) = 2a1S + 2a2z + 2ao E=2 Uj. Here, fi(') are linear functions, and fi(w) have the form: ki fi(t) = -E m,, I ki < n, i= 1,N. (10) m=li Define value I,(u) as the minimal value of 9p when constraints (8) hold. I),(-) is the confidence level a quantile function of random variable i's distribution, where ( = (I(u, w) = max{<)i(u, w)}. Then the initial problem can be reduced to the following form: ua = arg min < (u). (11) uEU 4. Guaranteeing solutior:gorithm. The quantile optimization problem (11) is computationally difficult, because the quantile function,a(u) can not be expressed analytically. Let us consider an algorithm for obtaining a guaranteeing solution of this problem. We use the General Minimax Approach ( see Malyshev and Kibzun 1987) to rewrite problem (11) in the equivalent form: (ua, E) = arg min max ((u, w), (12) uEU,EaEEa wEEa 6

where Ea is the class of confidence sets whith probability measure a. Denote B as the Borel v- field on the disturbance space. Then E, E B. Choosing an optimum set Eo E Ea is a difficult operation. However if we want to find a guaranteeing solution, we can replace Ea by some initial approximation Eo E En. For standard Gaussian random disturbances it is convenient to select Eo as a sphere Ca (see Malyshev and Kibzun 1987). The center of Ca is the mathematical expectation point of the random vector w and P(Ca) = a. In this case then the radius of Ca can be obtained analytically. Then the function:,(un) = max <(u, w) (13) wECa will be a majorant of the quantile function (a(u) (11), i.e.: a(u) < _a (u), Vu E U. (14) Given u E U and (13), oa(u) is the minimal value po, that satisfies all of the inequalities: max (i(u, w) <_ A, i =,N + 2. (15) wECa For every i = 1, N + 2, the maximum in (15) is attained at the point wt E ", which is a tangency point of the hyperplane fi(w) = const and sphere C, and can be found analytically. Hence, for calculating the guaranteeing solution, it is sufficient to solve the following linear programming problem: u = argmin o (16) uEU 7

subject to linear constraints: i,,i(uw?) _, i= 1,N (17) 'N+j(U) < c, j= 1,2. This problem can be solved by the standard simplex algorithm. Suppose that the guaranteeing solution (ius, Da) is discovered. Let us assume also that optimal area of solar powered plant Sac > 0 (the case, Sa = 0, will be discussed separately). Then the set of all vectors w satisfying the following inequalities: ' I(ia,w)~ _I, i1, -N, (18) is a polyhedral set: C = {w: fi(w) < ci,((, ), i=1,N}, (19) where - f(U~) -,Ci(a,) = ): i = 1,N, (20) a and Z = (i,w22,..., an). By construction of wi, the polyhedron fi(w) fi(w)), i = 1,N, (21) circumscribes sphere C,. By (21),(17) and that S > 0, all w which satisfy (21) belong to C. Hence Ca E C, and P(C) > a, because polyhedron (21) circumscribes the sphere Ca. If Sa = 0, then all feasible vectors w satisfy inequality (18), because in real problems probability distribution of random vector w has a finite support, and, so, the functions 8

fi(w), i = 1,N, have a finite value on the set of all feasible w. In this case, P(C)= 1. Let us consider a new confidence sphere C,1 with a radius Rai < R,, where Ra is the radius of C,. It is obvious that P(Cai) < a. Choose at the sphere Cal points w,* by analogy with the selection of w'. Now, consider problem (16) subject to: 4 (u, w*) <p, i= 1, N (22) ~N+j(U) < C, j 1,2. Let (As, <) be a solution to this problem. Denote: C- = {w: fi(w) < cic,(, Ua), i = 1,N}, (23) where the structure of ci(.) is described in (20). Lemma 1 If P(C*) > a, then Pa is an upper bound on the quantile function optimal value of problem (11) and the following inequality holds:,a < _a. Proof. 1. Consider set C* (23) as a confidence set insted of Ca. Then inequalities (17) are: Ad( U^~ S + MfiU~) < W i = 1,N, N(24) where u = (z, u2,...,u), and au = (z U,...,.... ~ = (~;, A*,.,). It is obvious, that (24) holds for the solution (Uti, a). Hence (Ui, Pa) is a feasible 9

solution for problem (16),(17), where w!, i = 1, N, are selected on the confidence set C*. It follows from the condition of the lemma, that P(C*) > a. As mentioned above, in the case when S, = 0, by analogy with the set C, P(C*) = 1, and (u., a<) is a feasiable solution, because of the finite value of the maximum fi(w), i = 1, N, on the set of all feasible w. Therefore, A, is an upper bound on the objective function optimal value in problem (11) (see Malyshev and Kibzun 1987). 2. Functions fi(') in (10) are linear combinations of Wj, j = 1,n. So by Ra > Rl1, the following inequalities hold: fs(wT*) < f,(wt), i =, N. (25) In fact, (ua, &a) is the solution of the problem (16),(17). So by (9),(25) and S > O0 inequalities (U wt*), i = 1, N + 2, (26) also hold. Hence, (&u, <ala) is a feasible solution for problem (16),(22), where the w,** are selected on the confidence set Cai. However (u., <) is the optimal solution of this problem. Therefore, Ia < 4O. This completes the proof of this lemma. O This lemma allows us to present the following algorithm for obtaining an improved upper bound on I,(.): Step 1. Find a solution (Ua, %,) of the linear programming problem (16),(17), where w', i = 1, N, are selected on the confidence sphere C, with probability measure a, and the w? maximize i(), i = 1, N, (9). Step 2. Check inequality P(C) > a, where C is from (19). 10

Step 3. If P(C) > a, then decrease the confidence sphere radius with the given increment ha and continue the cycle from Step 1. Step 4. In the opposite case, we select the previous step's linear program solution, bc,(), as an upper bound on the optimum value of ^a(') in problem (11) and set Ua = Uc. Remark. In the problem under consideration we approximate the true distribution of the random vactor w by a Gaussian distribution with given mathematical expectation and covariance matrix. All fi(w), i = i = 1, N, are linear functions, and so, for checking P(C) >_ a, it is sufficient to check the following inequality: F(c) > a, where c = (cl(>oa,ui ),..., cN(N a,u,)), and F(.) is a multidimensional Gaussian distribution function of the random vector (fi(w),..., fN(w)). For special methods for calculation F(c), see Deak 1980. In the case, So = 0, it is assumed, that P(C) = 1 and algorithm continues. 5. Numerical example and applying the quasi gradient algorithm for asymptotic exact solution. The considered problem is a system with n = 6 months. The results are obtained for two reliability levels, a = 0.99 and ac = 0.999. Initial data for this problem are: al = 3.75[unit/m2]; a2 = 10[unit/m3]; ao = 25[unit/m3]. The object water demand 11

per month: h = [29.6; 0.0001; 23.9; 36.2; 82.1; 173.4]. We approximate the distribution of the random vector w by a normal distribution with the following parameters: Mt = [0.00837; 0.00828; 0.0185; 0.0631; 0.123; 0.137] at = [0.000582; 0.000552; 0.00123; 0.00421; 0.00818; 0.00916]. Random variables wj, j = 1,6, are independent. We choose a short system life for obtaining compact results, that are convenient for analysis. The results are sensitive to the ao, al, a, Mt, at parameters relations. We omit the sensitivity analysis of this problem here. It is more important that this example reflects correlations between the guaranteeing and asymptotic exact solution of this problem, and the process of improving the guaranteeing solution using the guaranteeing solution algorithm. The equivalent problem is an unconstrained quantile optimization problem. This allows us to use the stochastic quasi gradint procedure for the quantile function minimization (see Kibzun and Kurbakovskiy 1991a, Kibzun and Kurbakovskiy 1991b) to obtain an asymptotic exact solution. This procedure is: uk+ = pu(uk- 0k(u)pk), u E U, (27) where pu(u) is projector to the set U; Pk is a nonrandom sequence of step multipliers; k is the number of the procedure step. The quantile function 4%i(u) stochastic quasi 12

gradient qk(u) is the following random vector: 1 n+1 0k(u) = 1 E[T((Ul, -- Uj + k,, * +l) - (Ul, -,j - Pk,., Un+l]ej, (28) where P3k is a nonrandom decreasing sequence term; n + 1 is the dimension of vector u; Uj, j = 1, n + 1, are random variables, which have uniform distribution on [uj - 13k, uj + 13k]; '(U) are independent quantile function estimates, obtained using a r sized sample of random variable I((u, w); ej, j =, n + 1, are basis vectors. To find the stochastic quasi gradient, qk(u), we propose to use the following statistics of the random variable 4((u, w): 4r,(u) = i[,ra] (29) c2,(u) = Ad - (,I - TIr-)(- + ln(T) + ln(1 - r)) where [Tra] is the whole path of the ra value; '[ra] is the [rca]th term in the variational sequence of the random variable D(u, w) sample, which have size r; (,Dr _r-i are extreme right terms of the variational sequence; p is the Euler constant. For the '1i(u) statistic, it is nesessary to use a sample size at least several times more than Ta = [ L-] +1, and to build the empirical cumulative distribution function of the random variable &(u, w). On the other hand for &r2(u), it is sufficient to use the sample size r = T, and we need only the extreme members of the sample to obtain this statistic. So 4r2(u) allows us to estimate high confidence level quantile function using small size samples. Statistical propeties of the estimates (28) were analysed in Kibzun 1991, Kibzun and Kurbakovskiy 1991b, Kurbakovskiy 1989. Conditions of the strong convergence 13

of uk to u, are considered in Kibzun 1991, Kibzun and Kurbakovskiy 1991a. Results of the guaranteeing and asymptotically exact solutions are given with the accuracy 0.01 for a = 0.99 and a, = 0.999 in Table 1. Table 1. Guaranteeing and Asymptotic Exact Solutions Guaranteeing solution Asymptotic exact solution ___ GS _ AES abs( GS-AES ) lD 4905.96 4995.02 4879.05 4973.78 26.91 21.24 S 1000.19 1017.03 1001.29 1017.22 1.1 0.19 V 61.33 63.21 59.14 61.11 2.19 2.1 w1 22.48 22.61 22.06 22.76 0.42 0.15 t2 0. 0. 0. 0. 0. 0. W3 0. 0. 0. 0. 0. 0. W4 O. 0. 0. 0. 0. 0. w O 0. 0. 0. 0. 0. 0. w6 O. 0. 0. 0. 0. 0. The process of improving the guaranteeing solution for a rithm from the section 4, is reflected in the Table 2. = 0.99, using the algo Table 2. Improving of the guaranteeing solution. r Guaranteeing solution P(C) 4.093 5189.61 0.9998 3.693 5100.92 0.9996 3.293 5015.09 0.9993 2.893 4931.99 0.9937 2.700 4905.96 0.9901 The probability of the set C was calculated using a Monte Carlo simulation procedure with sample size 10000. Befor the solving this problem, we transformed it to a problem, in which the random vector w has a standard normal distribution. So r in the Table 2 is the radius of the confidence sphere for a standard normal vector 14

(N(0, I)). We start the guaranteeing solution algorithm from the confidence sphere C, (a = 0.99) with radius r = 4.093. From Table 2, we see that in this case, using the guaranteeing solution algorithm, we can improve the standard guaranteeing solution (5189.61) by approximately 6 percent. Table 1 shows, that for the given conditions, the optimal water supply system satisfies demand by using fresh water produced by the solar powered plant. Only in the first month, when the plant productivity is small, do we need to use an external water supply. This fact can also be explained by the small capital cost of solar powered plant construction in this model example. Taple 1 illustrates, that, with increasing a, the guaranteeing solution tends to the asymptotically exact solution. This fact allows us effectively to use the guaranteeing solution algorithm by itself, and also to accelerate convergence of stochstic quasigradient methods using the guaranteeing solution as a good start point. 6. Conclusion. Choosing optimal water supply system structure is a typical optimization problem with probabilistic constraints. We presented nontraditional algorithm for this problem solution, reducing it to an unconstrained quantile optimization problem. After this reduction, the General Minimax Approach allows us to use linear programming to obtain an upper bound on the optimal value of the objective function, i.e the upper 15

bound on the cost for satisfying random water demand with given confidence level. An algorithm for improving this guaranteeing solution was also presented. This algorithm consists of a sequence of linear program solutions. Numerical tests showed, that an improved upper bound is close to the asymptotic exact solution, which was obtained using the stochastic quasi gradient algorithm for quantile function minimization. Moreover, numerical examples illustrate that the guaranteeing solution tends to an asymptotic exact one, when a increases. So the guaranteeing solution algorithm can be effectively used by itself and also for obtaining a good starting point for the stochastic quasi gradient algorithm. Combination of two considered methods is espesially useful for optimization in high reliability level systems. The approach outlined here could be extended to other similar optimization problem with probabilistic constraints. Acknowledgements - The authors are indebted to Dr. Vitali Yu.Kurbakovskiy and Dr. Bairam M.Muradov for many helpful discussions. References. Deak I.,Computation of multiple normal probabilities, in: Recents results in stochastic programming, ed. P.Kall and A.Prekopa, proceedings, Springer-Verlag (1980). Kibzun A.I.,Probabilistic optimization problem, Department of Industrial and Operations Engineering, Technical report 91-34, University of Michigan, Ann Arbor 1991. 16

and V.Yu.Kurbakovskiy, Numerical methods for runway space minimization under random disturbances,in Izv. AN USSR, Technic. Kibern., 1 (1991a) (in Russian) and V.Yu.Kurbakovskiy, Guaranteeing approach to solving quantile optimization problems, Anals of Oper. Res., v30 (1991b) Kurbakovsky V.Yu., Accelerate method for high confidence level quantile estimating, in collection of works VNIISI: Nonhomogeneous system dynamic, vol 14- VNIISI, Moscow,(1989) (in Russian) Malyshev V.V. and A.I.Kibzun, Analysis and synthesis of aircrafts high accuracy control, Mashinostroenie, Moscow, (1987), (in Russian) Mayer J., Probabilistic constrained programming: redused gradient algorithm implemented on PC, International institute for applied system analysis, A 23-61, Laxenburg, Austria,(1988). Mayer J., A nonlinear programming method for the solution of a stochastic programming model of A. Prekopa, Survey of Math. Progr., North Holland Publishing Co., New York, v 2 (1980) 129-139. Prekopa A., Numerical solution of probabilistic constrained programming problems, in: Numerical Techniques for Stochastic Optimization, ed. Yu.Ermoliev and R.J.B.Wets, Springer (1987). Network planning using two-stage programming under uncertainty, in: resent Resalts in Stochastic Programming (Proceedings of the international Contherence on Stochastic Programming, Germany,1979), Lecture Notes in Economics and Mathe 17

matical Systems 1979, Springer Verlag (1980),215-237. and T.Szantai, Floid control reservoir system design, Math. Progr. Stady, 9 138-151 and S.Ganezer, I.Deak and K.Patyi,The STABIL Stochastic Programming Model and its Experimental Application to the Electrical Energy Sector of Hungarian Economy, in: Stochastic programming, Proceeding of the International Conference of Stochastic Programming, Oxford, England,1974, (edited by M.A.N. Dempster), Academic Press, (1980),pp 369-385. and P.Kelle Reliability type inventory models based on stochastic programming,Math. Progr. Stady, 9 (1978),43-58. 18