Decision Support System Mine Problem Solver for Nonlinear Multi-Criteria Analysis

NOT FOR QUOTATION WITHOUT PERMISSION OF THE AUTHOR

DECISION SUPPORT SYSTEM MINE PROBLEM SOLVER FOR NONLINEAR MULTI-CRITERIA ANALYSIS

S. Kaden T. Kreglewski

January 1986 CP-86-5

CbLLaborative Papers r e p o r t work which has not been performed solely at t h e International Institute f o r Applied Systems Analysis and which has received only limited review. Views o r opinions expressed herein do not necessarily r e p r e s e n t those of t h e Insti- tute, its National Member Organizations, o r o t h e r organizations supporting t h e work.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS 2361 Laxenburg, Austria

PREFACE

The Regional Water Policies p r o j e c t of IIASA w a s focused o n inten- sively developed r e g i o n s w h e r e t h e water r e s o u r c e s are i n t e g r a t i n g ele- ments of t h e environment. The r e s e a r c h w a s d i r e c t e d t o w a r d s t h e develop- ment of methods a n d models t o s u p p o r t t h e r e s o l u t i o n of conflicts within s u c h socio-economic environmental systems. One of o u r case s t u d i e s d e a l s with open-pit lignite mining areas. The developed Decision S u p p o r t System MINE h a s b e e n implemented f o r a test r e g i o n in t h e Lusatian Lignite D i s t r i c t of t h e GDR.

The complex problems of s u c h r e g i o n a l policy a n a l y s i s are n o t tract- a b l e in o n e model using a n y of e x i s t i n g computational methods. T h a t is why a h e u r i s t i c two-level model a p p r o a c h h a s b e e n applied. Simplified first-level models t o g e t h e r with i n t e r a c t i v e p r o c e d u r e s f o r multi-criteria a n a l y s i s are used in t h e R u n n i n g Model f o r s c r e e n i n g analysis of r a t i o n a l long-term policies. Second-level models s e r v e f o r t h e v e r i f i c a t i o n a n d s p e c i f i c a t i o n of t h e r e s u l t s of s c r e e n i n g analysis.

In developing t h e system o u r major goal w a s t o make i t user-friendly, highly i n t e r a c t i v e a n d r o b u s t . F o r t h e planning model t h e s e f e a t u r e s are determined a b o v e a l l by t h e e f f e c t i v i t y of t h e problem s o l v e r f o r multi- c r i t e r i a analysis. The given p a p e r d e s c r i b e s s u c h a problem s o l v e r being developed f o r t h e DSS MINE. This r e s e a r c h h a s b e e n done within t h e frame- work of a c o l l a b o r a t i v e a g r e e m e n t between IIASA a n d t h e Technical Univer- s i t y of Warsaw, I n s t i t u t e of Automatic Control.

S e r g e i Orlovski P r o j e c t L e a d e r

Regional Water Policies P r o j e c t

The Decision S u p p o r t System MINE h a s been developed f o r t h e analysis of regional water policies in open-pit lignite mining areas. I t i s based on a two-level model a p p r o a c h . The first-level pLanning model is used f o r t h e estimation of r a t i o n a l s t r a t e g i e s of long-term development applying dynamic multi-criteria analysis. The second-level m a n a g e m e n t model considers managerial/ operational a s p e c t s f o r s h o r t e r time s t e p s (monthly and yearly).

The p a p e r d e s c r i b e s t h e problem s o l v e r f o r multi-criteria analysis in t h e planning model. This analysis i s based on t h e r e f i r e n c e p o i n t a p p r o a c h . For t h e solution of t h e resulting nonlinear programming problem t h e MSPN- algorithm, developed at t h e Institute of Automatic Control of t h e Technical University h a s been adopted. The s o l v e r considers t h e special c h a r a c t e r i s - t i c s of t h e mathematical model of t h e DSS MINE, as i t s non-linearity and t h e s p a r s e c h a r a c t e r of t h e resulting Jacobian matrix.

S t a r t i n g with t h e description of t h e g e n e r a l mathematical s t r u c t u r e of t h e planning model within t h e DSS MINE t h e problem formulation f o r multi- c r i t e r i a analysis based on t h e Refirence P o i n t A p p r o a c h i s given. Next, t h e non-linear problem s o l v e r MSPN i s p r e s e n t e d , including a program descrip- tion. Finally t h e r e s u l t s of some computational tests are shown.

CONTENTS

1. Introduction

2. Multi-Criteria Analysis

2.1 S t r u c t u r e of t h e planning model of t h e DSS MINE 2.2 Problem formulation f o r multi-criteria analysis 2.3 Scalarizing method

-

t h e R e f e r e n c e Point Approach 3. Non-linear Problem S o l v e r MSPN

3.1 T h e o r e t i c a l background

3.1.1 General d e s c r i p t i o n of t h e algorithm 3.1.2 Reduced g r a d i e n t algorithm

3.1.3 Penalty s h i f t algorithm 3.1.4 Verification of g r a d i e n t s 3.2 P r o g r a m d e s c r i p t i o n

3.2.1 P r o g r a m s t r u c t u r e

3.2.2 S t o r a g e method f o r Jacobian matrix 3.2.3 Optimization c o n t r o l p a r a m e t e r s 3.2.4 E r r o r handling

4. Computational Tests

4.1 Robustness of MSPN-solver

4.2 Influences of s t a r t i n g point values 4.3 Conclusions

R e f e r e n c e s Appendices

A COMMON blocks

B Subroutines and functions

DECISION SUPPORT SYSTEX

MINE

PROBLEM SOLVER FOR NONLKNEAR MULTI-CRITERIA ANALYSIS

S.

ad en'

and T. ~ r e g l e w s k i ~

1. I n t r o d u c t i o n

Regions with open-pit lignite mining are c h a r a c t e r i z e d by complex and s t r o n g i n t e r a c t i o n s in t h e socio-economic environmental system with s p e c i a l r e g a r d to water r e s o u r c e s . Caused by lignite mining, above a l l t h e n e c e s s a r y mine d r a i n a g e , o r i g i n a t e significant conflicts between d i f f e r e n t i n t e r e s t groups. For a detailed description of t h o s e problems see Kaden et al., 1985a.

Due t o t h e complexity of t h e socio-economic environmental p r o c e s s e s in min- ing areas, t h e design of regional water policies and water use technologies as well as mine d r a i n a g e c a n only b e done p r o p e r l y based on a p p r o p r i a t e mathematical models. From a c r i t i c a l analysis of t h e state-of-the-art of modeling in lignite mining areas i t h a s been concluded, t h a t above a l l methods and models are r e q u i r e d to s u p p o r t t h e analysis and implementation of r a t i o n a l l o n g - t e r m r e g i o n a l w a t e r p o l i c i e s in open-pit lignite mining a r e a s , t o a c h i e v e a p r o p e r balance between economic welfare and t h e state of t h e environment, Kaden et a l . 1985b.

Towards t h a t goal t h e r e s e a r c h of t h e Regional Water Policies p r o j e c t of IIASA, in collaboration with r e s e a r c h institutes in t h e GDR, and in Poland, in t h e period 1984-1985 was d i r e c t e d . One of i t s major p r o d u c t s i s t h e D e c i s i o n S u p p o r t S y s t e m MINE, see Kaden et a l . 1985a. Kaden 1986. The DSS MINE h a s been imple- mented f o r a test region in t h e Lusatian Lignite District in t h e GDR.

The analysis of regional water policies in mining r e g i o n s i s a problem of dynamic multi-criteria choice. An advanced system of decision a i d s i s needed which allows, Kadsn et a l . 1986:

')~nternational institute f o r Applied S y s t e m s Analysis Laxenburg, Austria

%wtitute of Automatic Control, Technical University of Warsaw, Poland

-

to consider t h e c o n t r o v e r s y among different water u s e r s and i n t e r e s t groups,

-

to include multiple c r i t e r i a some of which c a n not b e evaluated quantitatively,

-

_to t a k e into t h e account t h e uncertainty and t h e s t o c h a s t i c c h a r a c t e r of t h e system inputs as well as t h e limited possibilities to analyze a l l t h e decisive n a t u r a l and socio-economic p r o c e s s e s and impacts,

-

to o f f e r a set of decision a l t e r n a t i v e s , demonstrating t h e n e c e s s a r y trade-offs between d i f f e r e n t water u s e r s and i n t e r e s t groups.

A t p r e s e n t no mathematical methods are available or p r a c t i c a l applicable consid- e r i n g a l l t h e s e problems in one single model. Only time-discrete h i e r a r c h i c a l model systems c a n satisfy all requirements. Frequently a l r e a d y a two-level model h i e r a r c h y s a t i s f i e s most requirements. For t h e DSS MINE s u c h a two-level system h a s been realized.

The first-level model i s a Planning Model f o r t h e dynamic multi-criteria analysis f o r a relatively small number of p l a n n i n g p e r i o d s , j =1,

..

.,J as t h e time s t e p f o r principal management/technological decisions. Variable time s t e p s are used s t a r t i n g with o n e y e a r and increasing with time up to 15 y e a r s .

The planning model s e r v e s f o r t h e estimation of r a t i o n a l s t r a t e g i e s of long- t e r m systems development. These s t r a t e g i e s are s e l e c t e d by multi-criteria analysis.

The second-level Management Model i s applied f o r t h e simulation of systems behavior f o r a l a r g e r number of smaller management p e r i o d s (monthly and y e a r l y time steps). I t is used to analyze managerial decisions by t h e help of s t o c h a s t i c simulation and to v e r i f y r e s u l t s obtained with t h e planning model.

The DSS MINE is intended to b e highly interactively, user-friendly and robust.

The realization of t h e s e goals depends above all on t h e effectivity of t h e basic mathematical methods and models. One of t h e fundamental algorithms is t h e algo- rithm f o r non-linear multi-criteria analysis in t h e planning model.

The given p a p e r d e s c r i b e s t h e s o l v e r f o r non-linear multi-criteria analysis of t h e DSS MINE. I t h a s been developed in collaboration between IIASA and t h e Insti-

t u t e of Automatic Control of t h e Technical University Warsaw, Poland.

2. Multi-Criteria Analysis

2.1. Structure of the planning model of the DSS MINE

The planning model c o v e r s a p l a n n i n g h o r i z o n of 50 y e a r s divided into max- imum 10 planning periods, see Figure 1.

The f i g u r e i l l u s t r a t e s t h a t t h e highest a c c u r a c y is achieved f o r t h e f i r s t planning periods. The l a t e r planning p e r i o d s give rough estimates of f u t u r e systems development. T h e i r consideration e n s u r e s a r a t i o n a l systems development in t h e long-term run.

The planning model of t h e DSS MINE s e r v e s f o r t h e estimation of r a t i o n a l stra- t e g i e s of long-term systems development. These s t r a t e g i e s are s e l e c t e d by multi- c r i t e r i a analysis considering a number of c r i t e r i a . The c r i t e r i a have to b e chosen from a given set of i n d i c a t o r s , e.g. c o s t of water supply, cost of mine d r a i n a g e , satisfaction of water demand and environmental requirements. These indicators are assumed t o be i n t e g r a l values o v e r t h e whole planning horizon. In Figure 2 a block scheme of t h e planning model i s given.

With t h e purpose of a unified model being independent on t h e chosen c r i t e r i a i t i s assumed t h a t f o r all indicators bounds a r e given and all indicators are t r e a t e d as constraints. Based on t h a t t h e following multi-criteria problem for a s u b s e t OI l E L o of t h e indicators 0

(q

, l =I.

...,

L ) i s defined:

PLANNING HORIZON

4 b

SUBHORIZON

Figure 1: Time discretization f o r the planning model

Indicators of systems development O (

...,

l(i)~D(j),sv(j),sd(j),-*-)< m x 0

-... --..

L

Hydrological/socio- economic input l ( j )

Figure 2: Block schema of t h e planning model

-

b

j -1 j + 1

L

State of the system .Descriptive values

S,(j) = fS,(j* I* D.Sv)

Constraints C(j. Sv* S,, I)

<

0

v

b

State variables

Sv(j) = fSv(j,Sv(j

-

^l).D.Sd)

v v

-

^Decisions

D ( i ) ; 0,

0,

=

M i n i m u m ! I E L , s u b j e c t t o inequality c o n s t r a i n t s

C b ( j ) S O , j = l ,

...,

J equality c o n s t r a i n t s

bounds

This model d e s c r i b e s a non-linear dynamic multi-criteria problem. F o r t h e given problem in mining areas i t c a n b e assumed t h a t t h e systems dynamic is d e t e r m i n e d a b o v e a l l b y t h e e x t e r n a l l y fixed mine d r a i n a g e . The i n t e r n a l systems dynamic i s r e l a t i v e l y s l i g h t , t h a t means t h e influence of t h e state v a r i a b l e s on t h e r e s u l t ( i n d i c a t o r s ) i s l e s s important. Consequently t h e problem Eq.(2.1)-(2.4) may b e divided into s u b p r o b l e m s f o r a few s u b h o r i z o n s m , m =I,.

. .,

M, see F i g u r e 1.

o'(m)

=>

M i n i m u m ! , m =1,

...,

M (2.5) s u b j e c t t o Eq.(2.2)-(2.4) f o r s u b h o r i z o n m . This a p p r o a c h r e d u c e s t h e computa- t i o n a l e f f o r t d u e t o t h e s m a l l e r dimension of t h e non-linear programming problem.

In F i g u r e 3 t h e s t r u c t u r e of t h e Jacobian m a t r i x i s d e p i c t e d f o r a subhorizon with two planning p e r i o d s . The numbers give t h e a c t u a l s i z e of t h e problem f o r t h e

GDR

test area.

The F i g u r e i l l u s t r a t e s t h e s p a r s e c h a r a c t e r of t h e matrix. With t h e i n c r e a s i n g number of planning p e r i o d s p e r s u b h o r i z o n t h e matrix i s g e t t i n g more s p a r s e . The algorithm f o r non-linear programming h a s t o c o n s i d e r t h i s p r o p e r t y in o r d e r t o r e d u c e s t o r a g e consumption a n d computational e f f o r t .

2.2. P r o b l e m f o r m u l a t i o n f o r m u l t i - c r i t e r i a a n a l y s i s

I n s t e a d of t h e problem o r i e n t e d model formulation a b o v e f o r simplicity a n d convenience in t h e following a m o r e compact mathematical formulation of t h e multi-criteria problem Eq.(2.1)-(2.5) i s used.

W e c o n s i d e r a n o n l i n e a r optimization problems of t h e form:

minimize set of functions

s u b j e c t to:

nonlinear inequality c o n s t r a i n t s

nonlinear equality c o n s t r a i n t s

a n d bounds f o r a l l v a r i a b l e s

Figure 3: S t r u c t u r e of t h e Jacobian matrix

The nonlinear functions j i ( z ) , i =1,

...

^{, n h are assumed to} b e differentiable and t h e i r g r a d i e n t s must b e known in t h e analytical f o r m . These functions, t o g e t h e r with t h e r i g h t hand sides bi , i =1,

...,

^nh

,

t h e lower bounds l j

,

j =1,

...,

n

,

and t h e u p p e r bounds uj , j =l,.

. .

,n , constitute t h e model of systems behavior.

A s explained above some of t h e values j'{(z)

,

i 4 , .

. .

^, n, calculated in t h e model are defined as indicators of systems development. The set of indices I , con- tains t h e numbers of functions j'{(z) s e l e c t e d as objectives, being of i n t e r e s t f o r t h e decision maker applying t h e DSS MINE. This set c a n b e changed at any time.

Most frequently t h e set I , h a s more than o n e element and in such a case prob- lem (2.6)

-

(2.9) i s t h e problem of m u l t i o b j e c t i v e o p t i m i z a t i o n . To solve such a problem, o n e must s c a l a r i z e i t , i.e. r e d u c e i t to single c r i t e r i a equivalent using a s c a l a r i z i n g & n c t i o n .

2.3. Scalarizing method

-

the Reference Point Approach

For t h e DSS MINE t h e R e f e r e n c e Point Approach (Wierzbicki, 1983) is applied.

In t h i s method t h e reduction of t h e multiobjective optimization problem to a single objective one must b e interactively defined by t h e decision maker ( t h e model u s e r ) . The p r e f e r e n c e s among s e v e r a l c r i t e r i a are unknown a p r i o r i and are determined during t h e i n t e r a c t i v e ' p r o c e d u r e . For t h i s p u r p o s e , t h e decision maker defines a r e f e r e n c e value r{ f o r e a c h selected objective function j i ( z )

.

These values should r e f l e c t in some s e n s e d e s i r e d values of t h e objectives.

If all functions of t h e mathematical model are l i n e a r , then t h e scalarizing function using r e f e r e n c e values c a n b e given as follows:

with

w

=

^{(x )}

- ^-

distance between c u r r e n t values P i ( x ) and d e s i r e d ( r e f e r e n c e ) values ri ;

E

-

small positive number.

The second t e r m of t h e scalarizing function i s added t o g u a r a n t e e P a r e t o optimal- ity of t h e solution. (If E = O only weak P a r e t o optimality is guaranteed).

If P i ( x ) r r i f o r a l l objectives, t h e n t h e f i r s t t e r m of (2.10) is simply a Che- byshev norm of t h e distance between solution f ( ~ ) and r e f e r e n c e point r ; in a more g e n e r a l case, t h i s scalarizing function is not necessarily equivalent t o a norm, but i s s t r i c t l y monotonous and t h u s g u a r a n t e e s P a r e t o optimality of i t s mini- mas. Unfortunately, t h i s function is non-differentiable and cannot b e used in t h e case when t h e mathematical model is non-linear, because non-linear non- differentiable optimization methods are n e i t h e r r o b u s t n o r efficient enough f o r use in i n t e r a c t i v e systems. In t h e non-linear c a s e only a smooth approximation of function (2.10) c a n b e used.

The approximation used in t h i s package h a s t h e form:

with

no

-

number of objectives;

s i

-

scaling f a c t o r f o r i - t h objective;

W i

= Pi

(x )

-5 _- -

m e a s u r e of distance between

Pi

( x ) and ri ;

-

^Ti

-Pi

Pi -

lower bound f o r f i ( x ) and ri ;

P

-

positive (even) integer.

F o r t h e lower bound

~7

t h e u t o p i a (ideal) v a l u e may b e used minimizing objective i s e p a r a t e l y .

If p i s v e r y l a r g e , t h e n t h e approximation of t h e function (2.10) by (2.11) i s good. Unfortunately, t h e problem of minimizing (2.11) becomes badly conditioned in such a case. T h e r e f o r e , p = 4 o r 6 i s used in t h i s package.

3. Non-Linear Problem Solver MSPN 3.1. Theoretical background

3-1.1. General description of t h e algorithm

F o r a given fixed set of indices

I.

, fixed r e f e r e n c e point values ri and scal- ing f a c t o r s si , t h e resulting optimization problem t a k e s t h e form:

The feasible set X i s determined as a i n t e r s e c t i o n of t w o o t h e r s e t s :

X ,

i s t h e set d e s c r i b e d by nonlinear inequalities (2.7) and equalities (2.8).

XL

i s t h e set d e s c r i b e d by lower and u p p e r bounds (linear inequalities) (2.9).

Thus, t h e problem (3.1) is a s t a n d a r d nonlinear constrained optimization problem.

In t h i s package a double i t e r a t i v e p e n a l t y a l g o r i t h m is used f o r t h e problem solver.

The lower level algorithm solves t h e problem:

I 1

min F ( z ) = s , ( z ) + p ( z , v . k ) ,

D E xL

1

The objective function used h e r e , called penalty function, consists of the' sum of t h e original objective function and a penalty t e r m p ( z , v , k ) ( t h e p r e c i s e f o r m of t h i s t e r m is given l a t e r , see (3.6)). Linear c o n s t r a i n t s are satisfied at e a c h s t e p of t h e algorithm because a s p e c i a l method of reduced g r a d i e n t i s used f o r t h i s pur- pose. Nonlinear c o n s t r a i n t s , however, a r e violated and t h e penalty term in (3.5) is r e l a t e d t o t h i s violation.

The u p p e r level algorithm a d j u s t s p a r a m e t e r s v and k in t h e penalty func- tion t o satisfy nonlinear constraints. A t e a c h s t e p of this algorithm, t h e lower level problem (3.5) i s solved. However, t h e r e q u i r e d a c c u r a c y of i t s solution depends on t h e violation of nonlinear c o n s t r a i n t s . Nonlinear c o n s t r a i n t s are strongly violated in v e r y f i r s t i t e r a t i o n s of t h e u p p e r level algorithm and, t h e r e f o r e , t h e lower level problem c a n b e solved v e r y roughly.

3.1.2. Reduced gradient algorithm

The algorithm d e s c r i b e d h e r e and applied in t h e software p a c k a g e u s e s g r a - dient reduction, t h a t is, a n elimination of some g r a d i e n t components. If, at some point, t h e value of a p a r t i c u l a r v a r i a b l e zi is between i t s bounds I( and ui

,

t h e n this v a r i a b l e c a n b e e i t h e r i n c r e a s e d o r d e c r e a s e d . However, if this value i s equal to one of t h e bounds, s a y , to t h e u p p e r bound ui , t h e n t h i s v a r i a b l e c a n b e only d e c r e a s e d . In such a c a s e , negative values of t h e objective g r a d i e n t component c a n not b e a c c e p t e d and are set t o z e r o ; t h i s v a r i a b l e will remain unchanged in t h e n e x t direction of s e a r c h . This modified gradient is called r e d u c e d b e c a u s e some of i t s components are set to z e r o and i t acts only in s o m e s u b s p a c e of t h e s p a c e of a l l variables.

The algorithm begins by calculating t h e g r a d i e n t of t h e penalty function (3.5) at some s t a r t i n g point z 0

.

This g r a d i e n t is t h e n r e d u c e d in such a way t h a t a nonzero s t e p in t h e d i r e c t i o n of s e a r c h c a n b e performed inside t h e set

XL .

The step-size in t h i s d i r e c t i o n is calculated using q u a d r a t i c approximations in t h e line s e a r c h method. In t h e resulting point, t h e g r a d i e n t is calculated and r e d u c e d again. F i r s t d i r e c t i o n is just opposite t o t h e r e d u c e d g r a d i e n t (minus gradient), t h e n e x t d i r e c t i o n s are conjugate d i r e c t i o n s c o n s t r u c t e d on r e d u c e d gradients.

After s o m e number of i t e r a t i o n s t h e algorithm resets itself and uses minus g r a d i e n t direction again.

Following notation and symbols will b e used in t h e detailed description of t h e algorithm :

E

-

a c c u r a c y p a r a m e t e r given f r o m t h e u p p e r level algorithm;

k

-

i t e r a t i o n number;

m

-

number of conjugate direction.

The algorithm c a n b e c h a r a c t e r i z e d by t h e following s t e p s : lo Initialize: S e t k

=

0 and m

=

0

2' Calculate gradient: g k

=

VF(Z k ,

3O Gradient reduction f o r e a c h i =1,

. . .

, n :

if g!>0 a n d x f = ~ { o r g:<0 and x f = u ( then s e t g:=0.

If t h e resulting s u b s p a c e is d i f f e r e n t t h a n t h i s obtained in l a s t preceding iteration, set m = O

.

4' S t o p test: if ( I g k

I I

5 E t h e n s t o p

5' Calculation of new direction: if m = O o r m is g r e a t e r t h a n t h e number of nonzero elements in syk t h e n g e t a simple direction

dk

=

-syk

otherwise c a l c u l a t e conjugate direction using Polak-Ribiere algorithm:

dk

=

-gk + p d k - 1 with

6

' Direction check: if

<

d k , g k

>

^2: 0 t h e n set m

=

0 and go back t o s t e p 5' 7' Step-size limit:

TM

=

min max _'

d-4

_{- 4 k}

1

'&in

8' Line s e a r c h : find step-size ? such t h a t

/ ( x k + i d k )

=

min{(xk + r d k )

[0.ry

If i t fails ,i.e. ?

=

0 t h e n set m = O a n d go back t o s t e p 5' 9' Step:

x k + ' = x b

+

^{i d k , k}

=

k + I , m = m + 1 g o t o s t e p Z 0

The a c t u a l algorithm implemented in FORTRAN is much more complicated. I t includes many s a f e g u a r d s and i t t a k e s into account round off errors and finite a c c u r a c y of computations.

3.1.3. Penalty shift algorithm

The penalty t e r m in (3.5) h a s t h e following form:

with

k,

-

positive penalty coefficients, vi

-

penalty s h i f t s ,

v 4 non-negative f o r i = l ,

. . .

, n, , v{ unconstrained f o r i =nu + l , .

..

^, nh

.

S t a n d a r d penalty a l g o r i t h m s u s e a n y method of unconstrained optimization t o solve (3.5). The p e n a l t y c o e f f i c i e n t s are t h e n i n c r e a s e d a c c o r d i n g t o t h e violation of c o n s t r a i n t s o b t a i n e d a n d (3.5) i s solved again. This p r o c e d u r e i s r e p e a t e d until t h e solution of (3.5) is f o r c e d by t h e penalty t e r m to b e sufficiently c l o s e to t h e feasi- b l e set. Most f r e q u e n t l y , t h e penalty c o e f f i c i e n t s become v e r y l a r g e which makes problem (3.5) ill conditioned.

In t h e s h i g t e d p e n a L t y jgunction a t g o r i t h m t h e r e is no need t o i n c r e a s e p e n a l t y c o e f f i c i e n t s as s t r o n g l y as in s t a n d a r d penalty algorithms; p e n a l t y s h i f t s are used i n s t e a d to i n c r e a s e t h e p e n a l t y e f f e c t . The p e n a l t y function is s h i f t e d in t h e d i r e c t i o n o p p o s i t e t o t h e c o n s t r a i n t violation. In t h e case of inequality con- s t r a i n t s , t h i s l e a d s to a s h i f t "inside" t h e admissible set XN ; t o s h i f t a c o n s t r a i n t f { ( x ) 5 b { inside, o n e must d e c r e a s e t h e r i g h t hand s i d e b{ o r , equivalently, set a positive value of t h e s h i f t p a r a m e t e r v{ in (3.6). The p e n a l t y t e r m becomes t h e n a c t i v e in a band measured by v{ along t h e c o r r e s p o n d i n g boundary of t h e set

XN .

In t h e case of equality c o n s t r a i n t s , penalty s h i f t s c a n b e e i t h e r positive o r nega- tive; t h e a d e q u a t e s h i f t is just in t h e d i r e c t i o n o p p o s i t e to t h e c u r r e n t violation of c o n s t r a i n t s when solving (3.5). In both c a s e s , p e n a l t y s h i f t s i n c r e a s e t h e r e l a t e d p e n a l t y t e r m in (3.6). Additional s a f e g u a r d s are employed to avoid s t o p p i n g t h e algorithm inside t h e f e a s i b l e set with r e s p e c t t o t h e a c t i v e c o n s t r a i n t .

The algorithm is c h a r a c t e r i z e d by t h e following s t e p s : 1' Set initial k{

=

k p a n d v{

= o ,

i =I,.

. .

^{, n h}

2' Solve problem (3.5) using t h e r e d u c e d c o n j u g a t e g r a d i e n t algorithm; c a l c u l a t e maximal violation of c o n s t r a i n t s at t h e solution point:

with

gv

-

norm of violation of inequality c o n s t r a i n t s : g v

=

max m a x ( O , ( f { ( x ) - 9 ) )

(-1. .... _"0

gf

-

norm of inequality c o n s t r a i n t s f o r c e d inside t h e f e a s i b l e set:

h v

-

i s t h e norm of violation of equality c o n s t r a i n t s : h v = max I f { ( x ) - b { I

{=nu+'

...

nh

3' S t o p test: if q i s l e s s t h a n a given a c c u r a c y coefficient 7 t h e n s t o p . 4' If i t i s a f i r s t i t e r a t i o n t h e n set c = q a n d go t o s t e p 6'

.

5' If q

>

d t h e n set p

=

d a n d g o t o s t e p 9' 6' P e n a l t y s h i f t s :

v{

=

max ( 0 , vc

+

f { ( x )

-

^{b i ) , i =1,}

. . .

, n u 7' If q > c t h e n set p = c a n d g o t o s t e p 9'

9' F o r e a c h c o n s t r a i n t violated m o r e t h a n p , i n c r e a s e t h e penalty c o e f f i c i e n t s ki =2

*

and d e c r e a s e t h e p e n a l t y s h i f t s vi = 0 . 5 * v$

.

Go t o s t e p 8"

.

In most cases c h a n g e s of penalty s h i f t s in s t e p 6' are sufficient t o g e t c o n v e r - gence. However, if t h e rate of c o n v e r g e n c e in not s a t i s f a c t o r y , p e n a l t y coeffi- c i e n t s are a l s o c h a n g e d in s t e p 9"

.

3.1.4. Verification of gradients

Depending o n t h e number of time p e r i o d s t a k e n i n t o a c c o u n t in t h e model, com- p a r e Section 2.1, t h e number of functions f $ ( z ) in t h e g e n e r a l form (2.6)

-

^{(2.9) of}

t h e model c h a n g e s from a b o u t t h i r t y t o s e v e r a l hundreds. The number of n o n z e r o g r a d i e n t elements c h a n g e s f r o m s e v e r a l h u n d r e d s t o s e v e r a l thousands. The model of t h e system itself i s r a t h e r complicated: i t i s n o t a single FORTRAN s u b r o u t i n e , b u t r a t h e r l a r g e set of i n t e r c o n n e c t e d s u b r o u t i n e s a n d d a t a blocks, see Kaden 1986. Thus i t i s v e r y e a s y t o make a mistake calculating a n a l y t i c a l f o r m s of t h e d e r i v a t i v e s of complex e x p r e s s i o n s in t h e model.

Unfortunately, optimization a l g o r i t h m s are v e r y s e n s i t i v e t o s u c h mistakes a n d become inefficient o r e v e n f a i l if t h e c h a n g e s of o b j e c t i v e v a l u e s a n d con- s t r a i n t s are inconsistent with t h e i r g r a d i e n t s .

T h e r e f o r e i t i s n e c e s s a r y t o c h e c k t h e consistency of a l l g r a d i e n t s a f t e r e a c h modification of t h e model. F o r t h i s p u r p o s e , a s p e c i a l numerical algorithm was included in t h e optimization p a c k a g e . In t h i s algorithm g r a d i e n t s are c h e c k e d numerically by applying a f i n i t e d i f f e r e n c e method.

According t o t h e L a g r a n g e t h e o r e m in a single dimensional case t h e t e r m : f'(z2) - f ( z 1 >

x2 - 2 1

i s e q u a l t o t h e d e r i v a t i v e of t h e function f' ( z ) at some point between z , a n d r e . Typical a l g o r i t h m s of g r a d i e n t estimation assume t h a t t h i s t e r m i s t h e b e s t a p p r o x - imation of t h e g r a d i e n t in t h e c e n t e r of t h i s i n t e r v a l . However, s u c h a p p r o a c h r e q u i r e s 2

*

n p o i n t s to e s t i m a t e a g r a d i e n t in t h e n dimensional case

.

A r e g u l a r simplex method c a n b e applied t o minimize t h e number of points w h e r e t h e functions f $ ( z ) h a v e t o b e calculated. At a neighborhood of s u c h a point z o a r e g u l a r simplex of n

+

1 v e r t i c e s z,, z z , .

..

, z , + ~ i s c o n s t r u c t e d with x in t h e c e n t e r :

In a r e g u l a r simplex a l l d i s t a n c e s r

=

(

1

z -zj ⁽ , j =l,.

. .

^,n + l a r e e q u a l and value r i s called t h e r a d i u s of t h e simplex. The problem i s t o find t h e d i r e c t i o n s

d j = z j - x o t h a t s p a n t h e r e g u l a r simplex.

The algorithm of g r a d i e n t v e r i f i c a t i o n c o n s i s t of two p a r t s : c o n s t r u c t i o n of a simplex ( s t e p s 1"

-

4 " ) a n d calculation of g r a d i e n t e s t i m a t e ( s t e p s 5'

-

^{7' ).}

1' Calculate s c a l a r s :

2' Calculate f i r s t direction:

3' Calculate n e x t d i r e c t i o n s r e c u r s i v e l y f o r j

=

2,.

. .

, n +1

but f o r e a c h d additionally changing only component number j -1 d l - '

=

- ( j -1) d l - '

4' Calculate v e r t i c e s of t h e simplex as:

zj = z o + d j , j = l ,

. . . ,

n + l 5' A t e a c h point z j , j = l , .

. .

^, n + 1 , calculate values

PI =

Pc(zj)

6' F o r e a c h functions P i ( z ) calculate a n estimate of i t s g r a d i e n t component:

7' Compare t h e s e numerical estimates with t h e gradients calculated analytically at point z,.

3.2. Program description

3.2.1. Program structure

The problem s o l v e r MSPN f o r non-linear multi-criteria analysis i s embedded in t h e complex DSS MINE as i t i s illustrated in Figure 4. A detailed description of t h e model system i s given in Kaden, 1986.

The MSPN package i s a set of 2 1 interlinked FORTRAN subroutines and functions.

In Figure 5 t h e s t r u c t u r e of t h e MSPN p a c k a g e i s depicted.

A detailed description of t h e used COMMON blocks and of a l l subroutine and func- tions is given in t h e Appendices A and B.

The MSPN package contains only t h e algorithms f o r multiobjective optimiza- tion. All input d a t a and t h e mathematical model (Eq.(2.6)-(2.9)) are p r e p a r e d out- side ouf MSPN in t h e DSS MINE. The same holds t r u e f o r t h e s t o r a g e , processing and output of optimization r e s u l t s . A detailed description i s given in Kaden, 1986.

The only link between model u s e r s and t h e MSPN package t a k e s place in t h e case of modification of optimization c o n t r o l p a r a m e t e r s (default values are defined), and in case of numerical problems. In Section 3.2.3 and 3.2.4 some infor- mations are given.

I

^{DATA BASE} A

^I

l NTACT MINWAT GRAPH

Interactive Main program of the Display

data w Decision Support of colour

display Model System

editing

SYSTEM 3 SIMULA OPTlM

Second level : Simuletion Estimation of Estimation of

submodels of eff icient sol. utopia ml.

systems state Monte Carlo for selected for selected

development simulation mult.uiteria criteria

BALANC

-

REGIST MSPN

Monthly systems Registration Non-linear programming

balancr of stochastic

results

-

STATEM MODEL PSTAVA

Monthly system Preparation of Preparation of

Psr-a objective functions

-

system delcriptivet

and constraints state functions

OBJECT CONSTR TRANS

generator functions

functions functiom

Figure 4: S t r u c t u r e of t h e Decision S u p p o r t System

MINE

3.2.2. Storage method for Jacabian matrix

S i n c e t h e Jacobian matrix i s s p a r s e (compare Section 2.1) i t s columns ( g r a - d i e n t s of c o n s t r a i n t s ) a r e s t o r e d in a s p e c i a l way using a n i n d i r e c t indexing method. The g r a d i e n t s of c o n s t r a i n t s are not s t o r e d as n-dimensional v e c t o r s (sequences of n elements) b u t r a t h e r as s e q u e n c e s of elements known t o b e 'active' i.e. such elements which c a n h a v e n o n z e r o values. The remaining elements ('non a c t i v e ' ) must b e known to b e e q u a l to z e r o during t h e e n t i r e optimization p r o c e s s .

The a c t i v e elements a r e assumed t o b e listed as s e q u e n c e s of elements o r d e r e d a c c o r d i n g t o t h e i r p l a c e in t h e o r i g i n a l n dimensional v e c t o r . This method of listing i s not o b l i g a t o r y b u t t h e p a c k a g e i s more efficient if a c t i v e elements are o r d e r e d in s u c h a way.

The beginning indices and t h e lengths of t h e s e s e q u e n c e s are s t o r e d in t h e m a t r i x of indices icon(). The n o n z e r o g r a d i e n t elements themselves a r e s t o r e d in t h e Jacobian area of a g e n e r a l p u r p o s e s t o r a g e matrix r ( ) ( s e e d e s c r i p t i o n s of C O M M O N blocks /optc/ and /opti/ in Appendix A).

The a r r a y icon() i s used f o r a d d r e s s i n g t h e a r r a y icon() itself and t h e s t o r a g e a r r a y r ( ) . I t h a s t w o logical p a r t s : t h e f i r s t p a r t h a s ng

+

nh elements and contains f i r s t level indices f o r a d d r e s s i n g t h e second p a r t of icon(), t h e second p a r t con- t a i n s t h e d e s c r i p t i o n s of g r a d i e n t s of c o n s t r a i n t s . The length of t h e second p a r t of icon() is e q u a l t o ng

+

nh plus twice t h e number of s e p a r a t e s e q u e n c e s of elements

used f o r s t o r i n g a l l t h e elements of t h e c o n s t r a i n t g r a d i e n t s .

I

^JACDIM + PUTJAC

I

GO OPTDM

. -

Allocation of b Arithmetic

storage arms idrntif ication

GOOPT VERGRA

Figure 5: S t r u c t u r e of t h e MSPN p a c k a g e For a given c o n s t r a i n t number i :

icon( i )

-

contains t h e index of a n element in t h e second p a r t of icon(), where t h e g r a d i e n t description of t h i s p a r t i c u l a r c o n s t r a i n t begins, define ki

=

icon ( i );

icon( ki )

-

contains t h e minus index of t h e a r r a y r ( ) where t h e f i r s t a c t i v e element of t h i s g r a d i e n t is s t o r e d . All o t h e r s a c t i v e elements of t h i s c o n s t r a i n t a r e s t o r e d in t h e n e x t con- s e c u t i v e elements of t h e a r r a y r ( ) ;

icon( k i + l )

-

contains t h e number of t h e f i r s t a c t i v e element of t h e g r a d i e n t in t h e whole n dimensional v e c t o r ; if i t i s non- positive, t h e n t h i s c o n s t r a i n t h a s no a c t i v e elements at all (in case of dummy constraint) and i t i s t h e end of t h e description of t h i s constraint;

icon( k i + 2 )

-

contains t h e length of t h e sequence of a c t i v e elements (it c a n be equal t o one);

icon( k i + 3 )

-

contains t h e number of t h e beginning of t h e n e x t sequence of a c t i v e elements of t h e gradient; if i t i s positive t h e n element icon( ki+4 ) p e r f o r m s t h e s a m e r o l e as icon( k i + 2 ), otherwise i t i s t h e end of t h e description of t h i s con- s t r a i n t .

P a i r s number/length are r e p e a t e d as many times as r e q u i r e d t o d e s c r i b e a l l s e p a r a t e nonzero sequences of elements of t h e gradient.

3.2.3. Optimization control parameters

According t o t h e optimization algorithm d e s c r i b e d in Section 3.1 t h e following c o n t r o l p a r a m e t e r s a r e needed:

Range (rk):

Roughly estimated r a n g e of changes of v a r i a b l e s during t h e optimization pro- c e s s , scaling of v a r i a b l e s i s useful.

default: r k =l.

Norm (eps):

The s t o p test in t h e r e d u c e d conjugate g r a d i e n t algorithm c h e c k s whether t h e norm of g r a d i e n t s of t h e penalty function is l e s s t h e n t h e value of eps (denoted as E in t h e s t e p 4' of t h e algorithm in t h e Section 3.1.2), if eps is t o small s t o p with ip=4.

default: eps =0.1 Violation (eta):

The s t o p test in t h e penalty s h i f t algorithm c h e c k s whether all c o n s t r a i n t s are violated l e s s then 7 (denoted as 7 in t h e s t e p 3' of t h e algorithm in t h e Sec- tion 3.1.3)

default: e t a Penalty (penco):

The initial value of penalty coefficients in t h e penalty s h i f t algorithm (denoted as k: in t h e s t e p 1' of t h e penalty s h i f t algorithm

-

Section 3.1.3). A s a n esti- m a t e t h e r a t i o of g r a d i e n t s of t h e objective function t o g r a d i e n t s of con- s t r a i n t s should b e used.

default: penco =I.

Iterations (Ism):

Maximal number of i t e r a t i o n s (calculations of model values) default: Ism = I 0 0 0

rho:

P a r a m e t e r f o r scalarizing function of t h e r e f e r e n c e point method, see p in Section 2.3, Eq.(2.11)

default: r h o = 4

In Section 4.1 t h e influence of control p a r a m e t e r s on t h e optimization p r o c e d u r e is analyzed.

3.2.4. Error handling

The basic presumption f o r a successful optimization is t h e c o r r e c t model f o r - mulation, above all t h e analytical gradients. The , l a t t e r c a n b e c h e c k e d using t h e verification algorithm, see Section 3.1.4.

From t h e MSPN p a c k a g e t h e following i n t e r r u p t s are realized in case of possi- b l e e r r o r s :

-

Optimal solution not found b e c a u s e of t h e limit of i t e r a t i o n number (ip=3).

-

Optimal solution not found because of numerical errors

-

r e q u i r e d a c c u r a c y not attainable. In t h i s case t h e c o n t r o l p a r a m e t e r eps might b e increased. But, usually t h i s i n t e r r u p t indicates t h a t analytical formulas f o r functions o r t h e i r gradients are wrong (ip=4).

-

Optimization algorithms c a n not start because of too s m a l l s t o r a g e a r e a r e s e r v e d in t h e a r r a y r ( ) in t h e COMMON block /optc/ (ip=5). F o r s t o r a g e allocation see Kaden, 1986.

-

Optimization f a i l s because of t h e empty feasible set (ip=6). This might b e caused by t o s t r o n g c o n s t r a i n t s and bounds o r by model formulation e r r o r s . In case of t h e i n t e r r u p t s (ip=3,4,6) t h e model output c a n b e used t o localize possi- b l e e r r o r s . This output includes:

d e p s and eta

-

values compared on t h e s t o p test with e p s and eta parame- t e r s , respectively; d e p s i s t h e c u r r e n t norm of t h e penalty function, d e t a is t h e c u r r e n t violation of con- s t r a i n t s

ip

-

e r r o r condition p a r a m e t e r ; i p = 2 means "optimal solution found"

i t e r

-

t o t a l number of model calculations (calls f o r subroutine c r i c o n ) .

Lines with numbers at t h e beginning d e s c r i b e a c t i v e c o n s t r a i n t s with:

** -

t h e i r number (column in t h e Jacobian matrix), w

-

t h e c u r r e n t value of c o n s t r a i n t ,

v

-

t h e c u r r e n t penalty s h i f t , penal- t h e c u r r e n t coefficient.

The value I =(v + w ) * p e n a l is a n approximation of a Lagrange multiplier with r e s p e c t t o t h e scalarized objective function. It may b e used f o r a post-optimal analysis.

4. Computational T e s t s

4.1.

Robustness of

MSPN

solver

In o r d e r t o analyze t h e r o b u s t n e s s of t h e MSPN algorithm with r e s p e c t t o t h e optimization c o n t r o l p a r a m e t e r s , see Section 3.2.3, a s e r i e s of numerical t e s t s with t h e DSS MINE have been performed.

The tests have b e e n done f o r a planning horizon of 7 planning p e r i o d s . A s c r i - t e r i a t h e following had b e e n selected:

dev-m

-

Deviation municipal water demand/supply , dev-i

-

Deviation i n d u s t r i a l water demand/supply , cost-mi

-

Total mine d r a i n a g e c o s t ,

cost-m

-

Cost f o r municipal water supply, cost-i

-

Cost f o r i n d u s t r i a l water supply.

F o r e a c h c r i t e r i a t h e utopia point w a s s e l e c t e d as r e f e r e n c e point. A s t h e s t a r t i n g point o n e with significant deviation t o t h e e x p e c t e d solution w a s used in o r d e r t o r e a l i z e a l a r g e number of i t e r a t i o n s .

In F i g u r e 6 some r e s u l t s a r e d e p i c t e d illustrating t h e influence of c o n t r o l parame- ters on t h e r e s u l t s . Only t h o s e c r i t e r i a are shown which are s t r o n g l y a f f e c t e d by t h e p a r a m e t e r s . F o r t h e c r i t e r i a dev-m, dev-i, cost-m t h e influence i s almost negli- gible.

From t h e s e tests t h e following conclusions c a n b e drawn:

range (rk):

The influence of t h i s p a r a m e t e r on t h e numerical r e s u l t s i s negligible. The v a r i a - tions are l e s s t h e n 1%. But a t o small value a f f e c t s t h e number of i t e r a t i o n s signifi- cantly. According t o Figure 6 a values between 0 . 1 and 5 are r e a s o n a b l e .

violation (eta):

The influence of t h i s p a r a m e t e r i s a g a i n small, l e s s t h e n 1%. Smaller numbers of e t a i n c r e a s e t h e number of i t e r a t i o n s . A s a compromise eta=l0-' should b e chosen.

norm (eps):

A s e x p e c t e d t h i s p a r a m e t e r s t r o n g e r e f f e c t s numerical r e s u l t s a n d number of i t e r a t i o n s . Only values g r e a t e r / e q u a l 0.05 could b e chosen. F o r t h e value e p s =0.01 t h e r e q u i r e d a c c u r a c y h a s n o t b e e n a t t a i n a b l e . The r e s u l t s d e v i a t e in a r a n g e between maximum 5 and 10%

-

q u i t e a c c e p t a b l e from t h e p r a c t i c a l point of view.

F u r t h e r m o r e in e a c h case Pareto-optimality ( a t l e a s t locally) w a s a c h i e v e d , com- p a r e r e s u l t s f o r cost-mi and cost-i in Figure 6 c . A s a good compromise between a c c u r a c y and number of i t e r a t i o n s e p s =0.1 should b e chosen.

penalty (penco):

The initial penalty c o e f f i c i e n t h a s b e e n v a r i e d between 0.5 and 10. I t d o e s p r a c t i - cally n o t e f f e c t t h e numerical r e s u l t s . The influence on t h e number of i t e r a t i o n s i s small f o r values between 0.5 and 5. Only in t h e c a s e of p e n c o =10 t h e number of i t e r a t i o n s i n c r e a s e d . A s a good compromise p e n c o =I. i s p r o p o s e d .

Above t h e influence of optimization p a r a m e t e r s on t h e c r i t e r i a as i n t e g r a l parame- ters h a s b e e n analyzed. Another i n t e r e s t i n g question i s t h e i r influence o n t h e v a r i a b l e s (decisions).

The same p a r a m e t e r combinations as d e p i c t e d in Figure 6 h a v e b e e n analyzed with r e s p e c t to t h e i r impact on t h e v a r i a b l e s . Analogously to t h e c r i t e r i a t h e MSPN-algorithm i s a l s o v e r y r o b u s t with r e s p e c t to t h e v a r i a b l e s . The influence of varying range and violation i s almost negligible. The deviations are l e s s t h e n 5%, in most c a s e s e v e n l e s s t h e n 1% ( r e l a t e d to t h e value of t h e v a r i a b l e ) . Only in o n e case t h e deviations are s t r o n g e r . This i s d e p i c t e d in F i g u r e 7 f o r t h e v a r i a b l e s g ,,

a n d ad,,, as decisions on water allocation ( w a t e r quantity). F o r d e t a i l s o n t h e meaning of t h e v a r i a b l e s see Kaden et al., 1985a.

The solution f o r period 1 d i v e r s significantly in t h e case of range r k = 0 . 1 from t h e

Cost

-

^mi Cost

-

ⁱ

,

t t - - b

1370 -.-a 562

-

1 360

@

^560-

1340

Range

0.1 0.5 1.0 5.0 0.1 0.5 1.0 5.0

Norm = 0.1 ; violation =

a)

-.-.-,

^{Norm =} 0.5; violation =

0

Number of iterations

.- -... -.

⁴

...-.-

1 340

Violation Violation

1

o4

¹

o - ~

1

o - ~

¹

o - ~

1

o - ~

1 0-2 Range = 1 .O; norm = 0.1

0

Number of iterations b)

,.

^{Range = 1 .O; norm = 0.5 '}

*Iterations has been stopped

Cost

-

^mi _h ^Cost

-

ⁱ

Norm

Range = 1 .O; violation =

"Required accuracy not attainable

Figure 6: Influence of optimization c o n t r o l p a r a m e t e r s

Period 1 Period 4

vb,ex:

:,

q , b m ,

, ,

1

.o

Range Range

Norm = 0.5; violation =

Figure 7: Influence of optimization p a r a m e t e r s on v a r i a b l e s o t h e r solutions. For t h a t two r e a s o n s are seen:

-

t h e optimization is performed with r e s p e c t t o c r i t e r i a as i n t e g r a l values o v e r t h e whole planning horizon. Due t o t h e increasing time s t e p s of planning periods l a t e r planning p e r i o d s r e p r e s e n t s a l a r g e r p a r t of t h e c r i t e r i a , g e t a h i g h e r weight.

-

t h e p a r a m e t e r combination of a small norm esp=0.1 with a rough violation q=0.01 i s not v e r y reasonable.

Nevertheless t h e water balance is satisfied. The i n c r e a s e d %,,-value i s compen- s a t e d by a r e d u c e d qb ,,,

.

The effect of t h e norm is more significant as i t should b e e x p e c t e d from t h e r e s u l t s f o r c r i t e r i a , see Figure 6c. In Figure 8 t h e r e s u l t s f o r f o u r v a r i a b l e s are depicted. These are t h e v a r i a b l e s with t h e s t r o n g e s t deviations. For a l l o t h e r v a r i a b l e s t h e deviations are l e s s t h e n 5%.

The deviations are hardly t o b e explained. Probably t h e y are above a l l caused by f l a t objective functions with r e s p e c t t o t h e given variables. Small numerical devia- tions due t o d i f f e r e n t a c c u r a c y could lead t o d i f f e r e n t solutions.

4.2. Influence _of starting point values

The optimization problem t o b e solved is non-linear in most of i t s p a r t s . For such a complicated mathematical model as i t i s given f o r t h e DSS MINE i t is p r a c t i - cally impossible t o analyze analytically t h e p r o p e r t i e s of t h e objective function with r e s p e c t t o convexity and extremal values. The existence of local optima h a s t o b e expected. O r , with o t h e r words, t h e estimated solution i s not necessarily a n global optimal solution.

Two principle possibilities are available t o c h e c k t h e solution behavior with r e s p e c t t o local/global optima:

-

Application of a n optimization p r o c e d u r e resulting p e r definition in a solution being a global optimum. Such p r o p e r t y posses some random s e a r c h methods, e.g. ASTOP, Born 1985. The numerical e f f o r t of such methods i s extremely high, t h e i r applicability f o r t h e given problem is still under study.

I

I I I I I Year

1 2 3 4 5

*

^{6 10}

1 2 7 Period

1 I I I I I Year

1 3 4 5 6 10 7 Period

Range, = 1 .O; violation = norm = 0.05

-

0.1 ,,, 0.5

...

Figure 8: Influence of t h e norm o n v a r i a b l e s

-

Experimental analysis varying t h e s t a r t i n g points f o r optimization.

In t h e following a few r e s u l t s f o r a n experimental analysis are given.

The most logical way t o analyze t h e influence of s t a r t i n g point values i s t h e random selection of t h e s t a r t i n g points between u p p e r a n d lower bound. This h a s b e e n done using a random g e n e r a t o r f o r uniform d i s t r i b u t e d random numbers. F o r t h e tests t h e same c r i t e r i a as d e s c r i b e d in Section 4.1 h a v e been considered.

Results f o r s e l e c t e d c r i t e r i a are listed in Table 1.

The Table i l l u s t r a t e s t h a t only t h e c r i t e r i a cost -mi a n d cost -i significantly depend on s t a r t i n g point values. The maximal deviation i s in t h e r a n g e of 10%.

An i n t e r e s t i n g question i s , w h e t h e r t h e d i f f e r e n t s t a r t i n g points r e s u l t in dif- f e r e n t local optima, o r t h e r e s u l t s are simply d i f f e r e n t Pareto-optimal solutions. In Figure 9 t h e r e s u l t s f o r t w o c r i t e r i a are g r a p h i c a l l y i l l u s t r a t e d .

I t i s o u t of question t h a t only 11 tests are s t a t i s t i c a l l y not sufficiently f o r general- ization. N e v e r t h e l e s s from t h e Figure 9 could b e concluded t h a t some of t h e solu- tions (connected by t h e d a s h e d line) are global Pareto-optimal, a f e w o t h e r s only local. But e v e n f o r t h e "worst solution" t h e d i s t a n c e t o t h e n e x t h y p o t h e t i c global

Table 1: Solutions f o r randomly s e l e c t e d s t a r t i n g point v a l u e s

/ c r i t e r i a [ r e f

.-

utopia- nadir-1 solutions

7

[ I / sec.

1 )

d e v i ( O 1 1 0 1106

Cost

-

mi

1400

3 2 4 4 1 2 1 1 1 1 1 1

111 / sec.

1

cost --mi [ M i l l . M ] cost*

[ M i L L . M ] cost i [ M i l l . M ]

I I 1 b

530 540 550 560 Cost

-

ⁱ

F i g u r e 9: Solutions f o r d i f f e r e n t s t a r t i n g p o i n t v a l u e s 1151

12 509

Pareto-optimal p o i n t i s less t h e n 7

-

10% of i t s value.

Another s e r i e s of numerical tests h a s b e e n d o n e chosing lower a n d u p p e r bounds of v a r i a b l e s as s t a r t i n g p o i n t values. F o r t h e s e tests only o n e planning p e r i o d h a s b e e n analyzed. F o r e a c h v a r i a b l e o n e r u n was made with lower a n d u p p e r bound, t h e r e s u l t s were c o m p a r e d with a n a r b i t r a r y "nominal" solution.

F o u r t e e n v a r i a b l e s h a v e b e e n v a r i e d (28 r u n s ) . In Table 2 some r e s u l t s are dep- i c t e d .

The r e s u l t s i l l u s t r a t e t h e smal.1, influence of v a r i a t i o n s of single s t a r t i n g p o i n t s f o r o n e p e r i o d . The e f f e c t i s accumulating if m o r e planning p e r i o d s are u n d e r 1261 1261 1466 1315 1415 1411 1359 1443 1381 1446 1459 13 13 13 13 13 14 13 13 13 13 13 561 562 534 553 551 550 545 532 540 538 551 1151

12 509

1476 38 581

Table 2: S t a t i s t i c a l analysis of t h e influence of s t a r t i n g point on c r i t e r i a ( f i r s t p e r i o d )

1

c r i t e r i a

1

nominal value

!

mean value

1

s t a n d a r d deviation

1

consideration. In a l l cases f o r r e a s o n a b l e s t a r t i n g values t h e deviations are in t h e r a n g e of p r a c t i c a l a c c e p t a n c e .

d e v --m (1 / sec. ] d e v -i

[ L / sec. ] cost --mi

[Mill. Mark]

cost --m [Mill. Mark

1

cost i [Mill. Mark]

4.3. Conclusions

A s e r i e s of numerical tests h a s b e e n p e r f o r m e d in o r d e r t o analyze t h e r o b u s t n e s s of t h e MSPN-algorithm with r e s p e c t t o optimization p a r a m e t e r s and t o c h e c k t h e influence of s t a r t i n g point values. From t h e r e s u l t s c a n b e concluded t h a t t h e MSPN-algorithm i s well s u i t e d f o r t h e given problem. Variations in optimi- zation p a r a m e t e r s d o n o t e f f e c t t h e solutions significantly.

4 2 63.6

.

1.1 35.6

The mathematical model of t h e DSS-MINE b e h a v e s r o b u s t with r e s p e c t t o t h e selection of s t a r t i n g points. I t c a n not b e excluded t h a t local optima are estimated, b u t t h e local optima are e x p e c t e d c l o s e t o global optima.

F o r a l l tests deviations in c r i t e r i a values h a v e b e e n found l e s s t h e n a b o u t 10%.

Taking i n t o t h e a c c o u n t t h e a c c u r a c y of input d a t a , t h e simplified mathematical models of environmental a n d socio-economic p r o c e s s e s t h i s deviation i s fully a c c e p t a b l e .

4 2 64.2

1.1 35.5

Some of t h e tests are from t h e p r a c t i c a l point of view n o t r e a l i s t i c a l l y , espe- cially t h e random g e n e r a t i o n of s t a r t i n g point values. F o r t h e p r a c t i c a l problem t h e s t a r t i n g point values are in most cases known q u i t e w e l l a n d a consistency of v a r i a b l e s in time h a s t o b e considered. (It d o e s n o t make much s e n s e t o c h a n g e v a r i a b l e s as water allocation d r a s t i c a l l y between planning periods). Consequently t h e problem of s t a r t i n g point selection a n d local/global optima i s from t h e p r a c t i - c a l point of view less significant. F u r t h e r m o r e t h e r e s u l t s of t h e multi-criteria analysis within t h e p l a n n i n g model s e r v e s only as a guideline f o r a m o r e detailed analysis with t h e second level management model.

0.2 0.8 1.2 0.0 0.7

5. R e f ere nces

B0rn.J. 1985. Adaptively controlled random s e a r c h

-

a v a r i a n c e function a p p r o a c h . Systems analysis, Modelling, Simulation, Vol. 2, No. 2, pp.109-112, Akademie- Verlag Berlin.

Kaden,S., Hummel, J., Luckner.L., Peukert.D., Tiemer,K. 1985a. Water Policies:

regions with open-pit lignite mining (Introduction to t h e IIASA study), IIASA, WP-85-4, p.67.

Kaden,S., Luckner,L., Peukert,D., Tiemer,K. 1985b. Decision s u p p o r t model system f o r regional water policies in open-pit lignite mining areas. International Journal of Mine Water. Vol. 4, No.1, pp.1-16.

Kaden,S, Michels,I., Tiemer,K. 1986. Decision S u p p o r t System MINE; t h e Manage- ment Model, IIASA, CP-86, forthcoming.

Kaden,S. 1986. Decision S u p p o r t System MINE; Description of t h e Model System, IIASA, WP-86, forthcoming.

Wierzbicki,A.P. 1983. A Mathematical Basis f o r Satisficing Decision Making.

Mathematical Modeling USA 3:391-405 ( R e p o r t IIASA RR-83-7).