5th Conference on Optimization Techniques Part II [Part.2, 1 ed.] 3540066004, 9783540066002 [PDF]

Zitiervorschau

Lecture Notes in Computer Science Edited by G. Goos, Karlsruhe and J. Hartmanis, Ithaca Series: I.F.I.P.TC? Optimization Conferences

4 5th Conference on Optimization Techniques Part II

Edited by R. Conti and A. Ruberti

Springer-Verlag Berlin. Heidelberg New York 19 73

Editorial Board D. Gries • P. Brinch Hansen • C. Moler. G. Seegmtiller • N. Wirth Prof. Dr. R. Conti Istituto di Matematica "Ulisse Dini" Universit~t di Firenze Viale Morgagni 67/A 1-50134 Firenze/Italia Prof. Dr. Antonio Ruberti Istituto di Automatica Universith di Roma Via Eudossiana 18 1-00184 Roma/Italia

A M S Subject Classifications (1970): 6 8 A 2 0 , 6 8 A 5 5 , 9 0 A 1 5 , 9 0 B 1 0 , 9 0 B 2 0 , 92-XX, 9 3 A X X

I S B N 3"540-06600-4 Springer-Verlag Berlin • H e i d e l b e r g • N e w Y o r k I S B N 0-387-06600-4 Springer-Verlag N e w Y o r k • H e i d e l b e r g • Berlin

This work is subject to copyright. All rights are reserved, whether the whole or part of the material is concerned, specifically those of translation, reprinting, re-use of illustrations, broadcasting, reproduction by photocopying machine or similar means, and storage in data banks. Under § 54 of the German Copyright Law where copies are made for other than private use, a fee is payable to the publisher, the amount of the fee to be determined by agreement with the publisher. © by Springer-Verlag Berlin • Heidelberg 1973. Library of Congress Catalog Card Number 73-20818.Printed in Germany. Offsetprinting and bookbinding:Julius Beltz, Hemsbach/Bergstr.

PREFACE

T h e s e P r o c e e d i n g s a r e b a s e d on the p a p e r s p r e s e n t e d at the 5th I F I P C o n f e r e n c e on O p t i m i z a t i o n T e c h n i q u e s held in R o m e , May 7-11, 1973. The C o n f e r e n c e was s p o n s o r e d by the I F I P T e c h n i c a l C o m m i t t e e on O p t i m i z a t i o n (TC-7) and b y the C o n s i g l i o N a z i o n a l e d e l l e R i c e r c h e (Italian N a t i o n a l R e s e a r c h Council). T h e C o n f e r e n c e was devoted to r e c e n t a d v a n c e s in o p t i m i z a t i o n t e c h n i q u e s and t h e i r a p p l i c a t i o n to m o d e l l i n g , i d e n t i f i c a t i o n and c o n t r o l of l a r g e s y s t e m s . M a j o r e m p h a s i s of the C o n f e r e n c e was on the m o s t r e c e n t a p p l i c a t i o n a r e a s , including: E n v i r o n m e n t a l Systems, Socio-economic Systems, Biological Systems. A n i n t e r e s t i n g f e a t u r e of the C o n f e r e n c e was the p a r t i c i p a t i o n of s p e c i a l i s t s both i n c o n t r o l t h e o r y and in the field of a p p l i c a t i o n of s y s t e m s e n g i n e e r i n g . T h e P r o c e e d i n g s a r e d i v i d e d into two v o l u m e s . In the f i r s t a r e c o l l e c t e d t h e p a p e r s in which the m e t h o d o l o g i c a l a s p e c t s a r e e m p h a s i z e d ; in the s e c o n d those d e a l i n g with v a r i o u s a p p l i c a t i o n a r e a s . The I n t e r n a t i o n a l P r o g r a m C o m m i t t e e of the C o n f e r e n c e c o n s i s t e d of: R. Conti, A. R u b e r t i (Italy) C h a i r m e n , F e de Veubeke ( B e l g i u m ) , E. Goto (Japan), W. J. K a r p l u s (USA), J. L. L i o n s ( F r a n c e ) , G. M a r c h u k (USSR), C. Olech (Poland), L. S. Pontryagin" (USSR), E. R o f m a n ( A r g e n t i n a ) , J. S t o e r (FRG), J . H . W e s t c o t t (UK).

Previously published optimization conferences: Colloquium on Methods of Optimization, Held in Novosibirsk/USSR, (Lecture Notes in Mathematics, Vol. 112)

June 1968.

S y m p o s i u m on Optimization. Held in Nice, June 1969, (Lecture Notes in Mathematics, Vol, 132) Computing Methods in Optimization Problems. Held in San Remo, September 1968. (Lecture Notes in Operation Research and Mathematical Economics, Vol. 14)

TABLE OF CONTENTS

URBAN

AND

SOCIETY

SYSTEMS

S o m e A s p e c t s of U r b a n S y s t e m s of R e l e v a n c e to O p t i m i z a t i o n T e c h n i q u e s D; B a y t i s s

I

S e l e c t i o n of O p t i m a l I n d u s t r i a l C l u s t e r s f o r R e g i o n a l D e v e l o p m e n t S. C z a m a n s k i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

O p t i m a l I n v e s t m e n t P o l i c i e s in T r a n s p o r t a t i o n N e t w o r k s S. Giulianelli, A. La Bella .....................................

22

An

On-Line Optimization Procedure C. J. Macleod, A. J. A1-Khalili

for an Urban Traffic System .................................

Hierarchical Strategies for the On-Line Control of Urban Road Signals M. G. Singh ..................................................

51

Traffic

Application of Optimization Approach to the Problem of Land Use Plan Design K. C.: Sinha, A. J. Hartmann ....................................

42

60

Some Optimization Problems in the Analysis of Urban and Municipal Systems E. J. Beltrarni ~ ................................................ Combinatorial Optimization and Preference Pattern Aggregation J. M. Blin, A. B. Whinston ......................................

73

A Microsimulation Model of the Health Care System in the United States: The Role of the Physician Services Sector D. E. Yett, L. Drabek, M.D. Intriligator, L. J, Kimbell ............

85

C O M P U T E R AND C O M M U N I C A T I O N NETWORKS A Model for F i n i t e Storage M e s s a g e Switching Networks F. B o r g o n o v o , L. F r a t t a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

On C o n s t r a i n e d D i a m e t e r and M e d i u m O p t i m a l Spanning T r e e s F. Maffioli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

110

S i m u l a t i o n T e c h n i q u e s f o r the Study of M o d u l a t e d C o m m u n i c a t i o n Channels J. K. S k w i r z y n s k i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

118

~ p a p e r not r e c e i v e d

VI

Gestion Optimale d'un Ordinateur Multiprogramme ~ M6moire Virtuelle E. Gelenbe, D. P o t i e r , A. B r a n d w a j n , J . L e n f a n t . . . . . . . . . . . . . . . .

132

State-Space Approach in Problem-Solving Optimization A . S. V i n c e n t e l l i , M. S o m a l v i c o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

144

ENVIRONMENTAL

SYSTEMS

P e r t u r b a t i o n T h e o r y a n d t h e S t a t e m e n t of I n v e r s e P r o b l e m s G. I. M a r c h u k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

159

A M o d e l f o r t h e E v a l u a t i o n of A l t e r n a t i v e P o l i c i e s f o r A t m o s p h e r i c Pollutant Source Emissions R. A g u i l a r , L. F . G. d e C e v a l l o s , P . G. d e C o s , F . G 6 m e z - P a l l e t e , G. M a r t f n e z S a n c h e z . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . : .........

167

M a t h e m a t i c a l M o d e l l i n g of a N o r d i c H y d r o l o g i c a l S y s t e m , a n d t h e U s e of a S i m p l i f i e d R u n - O f f M o d e l i n t h e S t o c h a s t i c O p t i m a l C o n t r o l of a Hydroelectrical Power System M. F j e l d , S. L. M e y e r , S. A a m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

179

A Two-Dimensional M o d e l f o r t h e L a g o o n of V e n i c e C. C h i g n o l i , R. R a b a g l i a t i ....................................

203

Sea Level Prediction Models for Venice A. A r t e g i a n i , A. G i o m m o n i , A. G o l d m a n . n , P . S g u a z z e r o , A. T o m a s i n .................................................

213

Optimal Estuary Control Theory W. H u l l e t t

222

Aeration:

An Application

of D i s t r i b u t e d

Parameter

..................................................

Interactive Simulation Program for Water Flood Routing Systems F. Greco, L. Panattoni ..... : ............................

.....

An Automatic River Planning Operating System (ARPOS) E . M a r t i n o , B. S i m e o n e , T , T o f f o l i . . . . . . . . . . . . . . . . . . . . . . . . . . O n t h e O p t i m a l C o n t r o l o n a n I n f i n i t e P l a n n i n g H o r i z o n of C o n s u m p t i o n , Pollution, Population, and Natural Resource Use A. H a u r i e , M . P . P o l l s , P . Y a n s o u n i . . . . . . . . . . . . . . . . . . . . . . . . . . .

:.

251 241

251

ECONOMIC MODELS L i m i t e d R o l e of E n t r o p y i n I n f o r m a t i o n E c o n o m i c s J. Marschak ................................................. On a Dual Control Approach to the Pricing Policies Specialist M. A o k i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

264 of a T r a d i n g 272

VII

P r o b l e m s of O p t i m a l I n v e s t m e n t s w i t h F i n i t e L i f e t i m e C a p i t a l B. Nicoletti, L. Mariani ....................................... S o m e E c o n o m i c M o d e l s of M a r k e t s M. J . H. M o g r i d g e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . U t i l i z a t i o n of H e u r i s t i c s J. Christoffersen,

283 ...........

295

in Manufacturing Planning and Optimization P . F a l s t e r , E . S u o n s i v u , B. S v ~ i r d s o n . . . . . . . . . .

303

E c o n o m i c S i m u l a t i o n of a S m a l l C h e m i c a l P l a n t G. B u r g e s s , G. L . W e l l s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . An Optimal Growth Model for the Hungarian I. Ligeti

315

National Economy

324

.....................................................

BIOLOGICAL SYSTEMS O n O p t i m i z a t i o n of H e a l t h C a r e S y s t e m s J . H. M i l s u m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Theoretical and Operational Problems in Driving a Physical the Circulatory System B . A b b i a t i , R. F u m e r o , F . M . M o n t e v e c c h i , C. P a r r e l l a

335 M o d e l of ..........

347

M o d e l l i n g , S i m u l a t i o n , I d e n t i f i c a t i o n a n d O p t i m a l C o n t r o l of L a r g e Biochemical Systems J. P. Kernevez ..............................................

357

M o d @ l i s a t . i o n du T r a n s f e r t G a z e u x P u l m o n a i r e et C a l c u l A u t o m a t i q u e de la Capacit~ de Diffusion D. S i l v i e , H. R o b i n , C. B o u l e n g u e z . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

366

A Model for the Generation Movements A. Cerutti, On

of t h e M u s c l e F o r c e

F. Peterlongo,

R. Schmid

Some Models of the Muscle Spindle C. Badi, G. Borgonovo, L. Divieti

paper not received

~

During Saccadic Eye

.........................

.............................

_378

Contents

of

(Lecture

Notes

Part

I in

Computer

Science,

Vol.

3)

SYSTEM MODELLING AND IDENTIFICATION I d e n t i f i c a t i o n of S y s t e m s A. V. B a l a k r i s h n a n

Subject to Random State Disturbance ...........................................

1

Adaptive Compartimental Structures in Biology and Society R. R. M o h l e r , W . D . S m i t h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

O n t h e O p t i m a l S i z e of S y s t e m M o d e l M. Z . D a j a n i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

Information-theoretic Methods for Modelling and Analysing Large Systems R. E . R i n k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

A New Criterion for Modelling Systems L . W. T a y l o r , J r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

S t o c h a s t i c E x t e n s i o n a n d F u n c t i o n a l R e s t r i c t i o n s of I l l - p o s e d Estimation Problems E. Mosca ...................................................

57

Regression Operator Application to Some A. Szymanski An On

in Infinite Dimensional Vector Spaces and Its Identification Problems ............................................

Approach to Identification and Optimization in Quality Control W. Runggaldier, G. R. Jacur ................................... Optimal Estimation Y. A. Rosanov ~

Identification J. Cea The

69 83

Processes

de Domaines ......................................................

Modelling of Edible Oil Fat Mixtures J. O. Gray, J.A. Ainsley ......................................

DISTRIBUTED Free

and Innovation

...

92 103

SYSTEMS

Boundary Problems and Impulse Control J . L. L i o n s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A Convex Programming Method in Hilbert Space and Its Applications t o O p t i m a l C o n t r o l of S y s t e m s D e s c r i b e d b y P a r a b o l i c E q u a t i o n s K. M a l a n o w s k i . . . . . . . . . . . . . . . . : ..............................

~paper not received

116

124

About S o m e F r e e B o u n d a r y P r o b l e m s C o n n e c t e d with H y d r a u l i c s C. B a i o c c h i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

137

M~thode de D ~ c o m p o s i t i o n A p p l i q u e e au C o n t r S l e O p t i m a l de Syst~mes Distribu~s A. B e n s o u s s a n , R. G l o w i n s k i , J. L. L i o n s . . . . . . . . . . . . . . . . . . . . . . .

141

A p p r o x i m a t i o n of O p t i m a l C o n t r o l P r o b l e m s of S y s t e m s D e s c r i b e d by B o u n d a r y - v a l u e M i x e d P r o b l e m s of D i r i c h l e t - N e u m a n n T y p e P. C o l l i F r a n z o n e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

152

C o n t r o l of P a r a b o l i c S y s t e m s with B o u n d a r y C o n d i t i o n s I n v o l v i n g T i m e - D e l a y s P. K. C. Wang

..................................................

153

GAME T H E O R Y C h a r a c t e r i z a t i o n of C o n e s of F u n c t i o n s I s o m o r p h i c to C o n e s of Convex Functions J. - P . A u b i n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

174

N e c e s s a r y C o n d i t i o n s and S u f f i c i e n t C o n d i t i o n s f o r P a r e t o O p t i m a l i t y in a M u l t i c r i t e r i o n P e r t u r b e d S y s t e m J. - L . Goffin, A. H a u r i e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

184

A U n i f i e d T h e o r y of D e t e r m i n i s t i c T w o - P l a y e r s Z e r o - S u m Differential Games C. Marchal About

...................................................

Optimality of Time of Pursuit M. S. Nikol'skii ...............................................

194

202

PATTERN RECOGNITION A l g e b r a i c A u t o m a t a and O p t i m a l S o l u t i o n s in P a t t e r n R e c o g n i t i o n E. A s t e s i a n o , G. C o s t a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

205

A New F e a t u r e S e l e c t i o n P r o c e d u r e f o r P a t t e r n R e c o g n i t i o n B a s e d on Supervised Learning J. K i t t l e r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

218

On R e c o g n i t i o n of High D e f o r m e d P a t t e r n s by M e a n s of M u l t i l e v e l Descriptions S. T y s z k o * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A Classification Problem'in Medical Radioscintigraphy G. W a l c h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

250

T h e D y n a m i c C l u s t e r s Method and O p t i m i z a t i o n in N o n - H i e r a r c h i c a l Clustering E. D i d a y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

241

~ p a p e r not r e c e i v e d

XI

OPTIMAL

CONTROL

A Maximum Principle for General Constrained Optimal Control Problems - An Epsilon Technique Approach J . W. M e r s k y ................................................

259

O p t i m a l C o n t r o l of S y s t e m s G o v e r n e d b y V a r i a t i o n a l J. P. Yvon ...................................................

265

Inequalities

O n D e t e r m i n i n g t h e S u b m a n i f o l d s of S t a t e S p a c e W h e r e t h e O p t i m a l ValueSurface Has an Infinite Derivative H. L. S t a l f o r d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

276

C o n t r o l of A f f i n e S y s t e m s w i t h M e m o r y M. C. D e l f o u r , S . K . M i t t e r . . . . . . . . . . . . . . . . . . . . . .

292

...............

Computational Methods in Hilbert Space for Optimal Control Problems with Delays A. P . W i e r z b i c k i , A. H a t k o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

304

S u f f i c i e n t C o n d i t i o n s of O p t i m a l i t y f o r C o n t i n g e n t E q u a t i o n s V. I. B l a g o d a t s k i h .............................................

319

V a r i a t i o n a l A p p r o x i m a t i o n s of S o m e O p t i m a l C o n t r o l P r o b l e m s T. Zolezzi ...................................................

329

Norm P e r t u r b a t i o n

of

Supremum Problems

J. Barange r ................................................ On Two Conjectures P. Brunovsk~

333

about the Closed-Loop Time-Optimal ................................................

C o u p l i n g of S t a t e V a r i a b l e s i n t h e O p t i m a l L o w T h r u s t Transfer Problem R. H e n r i o n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . O p t i m i z a t i o n of t h e A m m o n i a O x i d a t i o n P r o c e s s of N i t r i c A c i d P . U r o n e n , E. K i u k a a n n i e m i .

.

.

.

.

.

.

.

.

.

.

Control

341

Orbital

345

Used in the Manufacture .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

360

STOCHASTIC CONTROL Stochastic Control with at Most Denumerable Number J. Zabczyk ..................................................

of C o r r e c t i o n s

370

D e s i g n of O p t i m a l I n c o m p l e t e S t a t e F e e d b a c k C o n t r o l l e r s f o r L a r g e Linear Constant Systems W . J . N a e i j e , P . V a l k , O. H. B o s g r a ............................

375

C o n t r o l of a N o n L i n e a r S t o c h a s t i c B o u n d a r y V a l u e P r o b l e m J . P . K e r n e v e z , J. P . Q u a d r a t , M. V i o t . . . . . . . . . . . . . . . . . . . . . . . . .

389

An Algorithm to Estimate Sub-Optimal Present Values for Unichain Markov Processes with Alternative Reward Structures S. D a s G u p t a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

399

XIi MATHEMATICAL PROGRAMMING S o m e R e c e n t D e v e l o p m e n t s in N o n l i n e a r P r o g r a m m i n g G. Z o u t e n d i j k . . . . . . . . . . . . ....................................

407

P e n a l t y M e t h o d s and A u g m e n t e d L a g r a n g i a n s in N o n l i n e a r P r o g r a m m i n g R. T. R o c k a f e l l a r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

418

On I N F - C o m p a c t M a t h e m a t i c a l P r o g r a m s R . J. -B. W e t s

426

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

~

.

.

.

.

.

.

.

.

, .

.

.

.

.

.

.

.

.

.

.

.

.

.

°

.

.

.

.

Nonconvex Quadratic P r o g r a m s , Linear C o m p l e m e n t a r i t y P r o b l e m s , Integer Linear Programs F . G i a n n e s s i , E. T o m a s i n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

. . o

and

457

A W i d e l y C o n v e r g e n t M i n i m i z a t i o n A l g o r i t h m with Q u a d r a t i c T e r m i n a t i o n Property G. T r e c c a n i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

450

A H e u r i s t i c A p p r o a c h to C o m b i n a t o r i a l O p t i m i z a t i o n P r o b l e m s E. Biondi, P. C. P a l e r m o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

460

A New Solution f o r the G e n e r a l Set C o v e r i n g P r o b l e m L. B. K o v ~ c s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

471

A T h e o r e t i c a l P r e d i c t i o n of the Input-Output T a b l e E. K l a f s z k y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

484

An I m p r o v e d A l g o r i t h m for Pseudo-~Boolean P r o g r a m m i n g S . W a l u k i e w i c z , L. S~omi~ski, M. F a n e r . . . . . . . . . . . . .

493

. . . . . . . . . . . .

N u m e r i c a l A l g o r i t h m s for Global E x t r e m u m Search J. E v t u s h e n k o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

~...

505

N U M E R I C A L METHODS G e n e r a l i z e d S e q u e n t i a l G r a d i e n t - R e s t o r a t i o n A l g o r i t h with A p p l i c a t i o n s to P r o b l e m s with Bounded C o n t r o l , Bounded State, and Bounded T i m e R a t e of the S t a t e A. M i e l e , J . N . D a m o u l a k i s ~ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G r a d i e n t T e c h n i q u e s f o r C o m p u t a t i o n of S t a t i o n a r y P o i n t s E. K. B l u m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

509

P a r a m e t e r i z a t i o n and G r a p h i c A i d in G r a d i e n t M e t h o d s J. - P . P e l t i e r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i .................

517

L e s A l g o r i t h m e s de C o o r d i n a t i o n dans la M~thode M i x t e d ' O p t i m i s a t i o n & Deux N i v e a u x G. G r a t e l o u p , A. T i t l i , T. L e f ~ v r e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

526

* p a p e r not r e c e i v e d

X!Ii

A p p l i c a t i o n s of D e c o m p o s i t i o n and M u l t i - L e v e l T e c h n i q u e s to the O p t i m i z a t i o n of D i s t r i b u t e d P a r a m e t e r S y s t e m s Ph. C a m b o n , L . Le L e t t y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

538

A t t e m p t to S o l v e a C o m b i n a t o r i a l P r o b l e m in the C o n t i n u u m by a Method of E x t e n s i o n - R e d u c t i o n E. S p e d i c a t o , G. T a g l i a b u e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

554

SOME ASPECTS OF URBAN SYSTEMS OF Pkw,~VANCE TO 0PTIMISATION TECHNIQUES - D. BAYLISS GREATER LO~DON COUNCIL

1.

Introduction

I would like to start by making clear what will shortly be obvious anyway - that I am not an expert in information science or optimisation techniques.

I come here

as an urban plsrmer to try and explain the nature of some of the problems to be faced in practice, in a way ~hich I hope is relevant to this important branch of technolo~ and also to try and learn something of the possibilities for the application of opti~isaticn techniques in tackling the difficult problems of planning and managing towns and cities. My theme will be that the essential characteristics of urban systems are so complex that they raise conceptual measurement and computational problems of an order which preclude the l~ssibility of building formal comprehensive optimisation models ~hich are of value to the practising planner.

There i~, however, a role for sectoral

optimisation models operated outside, but better, within a comprehensive analytical framework. 2.

Txpes of System Parameters

As I have already said urban planning is distinguished from many other types of planning by its complexity but also it is characteristically confronted with large arrays of actual and potential conflicts and the resulting judicial role presents one of the most difficult challenges for formal optimisation techniques. FIGURE 1 Polic,y Parameters in Urban Systems Physical

Locat ion

Social

Timing

Economic

Agencies?

Political? Structural change Management regimes Figure 1 shows some of the main parameters of interest to the urban policy maker. The physical, social and economic aspects of environmental change have to be considered, the way in which the structure of the city operates and might be developed under alternative management regimes, not simply, but differentiated in historic time and space. 0ptimisation models typically operate within one of these sectors, some of the more adventurous linking two areas but rarely all three.

Where the model deals with

structural and maaagement aspects or the dynamics of the system again they are

usually confined to a single sector. Another major problem is that of defining optlm~l conditions.

Urb~u systems are

typified by a multiplicity of client groups with conflicting and often undefined preferences.

Thus simple objective functions are rarely meauingful except in

single sector models with single client groups. FIGURE 2 Objective Parameters in Urbau Systems Feasibility Efficiency Equity Quality Also it is rarely possible or proper to predetemmine the weights to be attached to the criteria under these four heads (at least those which cau be measured) or to the interests of different client groups as those weights are traditionally determined by political rather than technical processes. 3.

Am Analytical Framework

Having suggested that there are both fundamental aud practical reasons why global optimisation models caunot be usefully developed for u~ban system the question must be faced as to how partial optimisation models should be most appropriately developed.

A major danger in the use of partial optimisation models is that they

attract too much weight to that part of the system with ~hich they deal and therefore, unless their scope dominates the system, their exercise can be counter productive to ~ood decision making.

This danger can be lessended by establishing

a comprehensive view of the stractu~e of the system as a context for identifying the scope of relevance of sectoral optimisation. Figure 3 outlines in simple diagrammatic form am analytical framework for an urban system.

Paradigms of this kind can reveal the partiality of particular models

and show linkages between individual sub-model.

The value of basic fr~eworks of

this kind can be increased by the development of behaviotu~al models of the relationships they portray and a number of such models exist.

Perhaps the best

known is that for the transportation system - illustrated in Figures 4 aud 5-

FIGURE 5 Relationship of Commercial Industrial and Recreational Centres

Production Employment

Consumption Retailing ~Education, etc.

Distribution ~ of wholesale_go_ods_.

/ Communications Infrastructure

..

',,% ~o//

• °

\~.

:|of ~us~

"to

!t service Centres

....

%

.

I/

~ !

~,

~.

imtionship of ,,'[i

I Housin~d

I W°rk Places .

• . . . . .

•

. . . . . .

. . . .

~

°o

Residence Housing

.....Spatial structure ~/So cio/economic structure --- Transportation stz~cture ~__-~External linkages

° . . .

FIGURE % I Infrastructure I

Spatial structure

,

1

Socio/economic structure

1 GENERATION

1 DISTRIBUTION

1 CHOICE

1 ASSIGNMENT

FIGURE 5 TRIP GENERATION Oi=

a+bx

i + cy i + azi

TRIP DISTRIBUTION Tij = OiDjAiGje -#Cij MODAL CHOICE -cij k Tij k = E E _cij* ROUTE CHOICE ~j =~nCij

T = 0 = D = i,j= C = k = R = x,y,z = a,b,c,d=

trips origins destinations zones cost mode route socio/economic variables coefficients

Models of this kind are often well developed and tested against extensive empirical evidence.

In addition to specifying important interactions between other sub-

models urban trs~sport models can provide the basis for sectoral optimisation models. There are many examples of this but one of the most interesting was developed in the London Transportation Study in order to devise an optimum traffic restraint policy - Figure 6. FIGURE 6 a b c C ~ i A'a + c~ I A'b + o< I A'c +

_Lc I

"+'oL n A'n

~Cx

x

There being 'n' zones and 'x' links n Subject to

i A'

% Ai i A'i ~-~-kAi

(Equity constraint )

A' i = attractions allowed in zone 'i' Ai

= attractions desired in zone 'i s

C. J

= capacity of link j

K

= constant

c~ i J

number of trips generated on link 'j ~ by a = unit attraction at zone

More extensive models exist ~ i c h crudely embrace aspects of the transport, spatial and socio-economic characteristics of Urban Systems.

Perhaps the best known is

that developed by Lowry originally for the Pittsburg area.

As can be seen from

Figures 7 and 8 this is a combined spatial structure/transportation structure built around the journey to work/service journey elements of the transportation system.

Simple optimisation can be based upon m~Ytmising choice or accessibility

or minimising transport costs.

,$ I Basic I Emoloyment IService Em~lo~yment I--

l LOWRY MODEL - STRUCTURE

FIG~m 7 ~

~ P°pulati°n dependant I on basic employment ~ Population generated I by service employment|

J

FIGURE 8 Tij = Ai(Eb) iPjae

~ij

-ucij

= (Eb)i

(Ebp) j = ~ i j Where Tij

= Basic work trips between zones i and j

Cij

= Cost of travel between zones i and j

(Eb)i = No. of basic jobs in zone i (Ebp)j= No. of basic workers in zone j a, U

= constant (calibrated)

Ai

=

Pj

= 'population' of zone j

I__~.Ipj ~"a e -ucij

LOWRY MODEL 4.

-

DISTRIBUTION OF EMPLOYEES IN BASIC ~MPLOYMENT

An Alternative Apprpach

What I would like to do now is to go on to illustrate a slightly different approach which allows an interactive process between formal optimisaticn and informal eptimisation ~hich is more congruent with the nature of the urban decision making process.

Thus non-controversial cptimisation can be buried snd treated exactly

within the model and controversial optimisation is exposed and treated subjectively outside the model. The model referred to is the General Urban Pl~a~uing model currently being developed by Dr. Young in the Department of Pls2Luing and Transportation of the Greater London Council. The number of variables and constraints makes it necessary to deal with the urban area by dividing into nine sections for which individual sub-models are operated a tenth sub-model reconciles the results of the sectoral sub-models. 5.

The Sector Sub-models

The model deals in terms of land uses (see Figure 9) aad their states through four consecutive time periods.

Each land use is described by four types of variables: Quantity Cost Population served Life

Also constraints of three types are attached to each land use excluding the simple logical constraints necessary to programming models; Compulsory

Arbitrary

these are

Balancing

FIGURE 9 LAND USES

(a) housing (b) hospit~s (c)

schools

(d) factories (e)

offices

(f)

shops

(g)

public open space

(h)

transport

(i)

vacant

(j) others 6•

Objectives

Objectives are set for housing density

cost, land and resource utilisation, etc.,

and the model run to maximise these subject to the constraints, the nature and value of the arbitrary constraints and the objective function are varied to expose conflicts, illuminate the scales of trade between objective elements as well as discerning composite optimum conditions using a linear programme. Although limited by data availability and the simplicity of assumptions employed this model is yielding valuable insights into the relationships between sectoral standards and objectives.

What is more the treatment of variables in both

formal end subjective modes creates a richness and sympathy ~hich minimises the probability of alienating the practicing u~ban planner. 7•

Conclusions

I am afraid that my comments have been too brief and patchy to do full justice to the nature of the problems of urban systems.

However, I hope they will help

a little to cl~rify my closing summary comments:(1)

The problems of optimisation in urban systems are created by:(a)

The need, in real life, to treat a rauge of complex sub-systems with quite different theoretical bases at one and the same time.

(b)

The inherent lack of concensus between client groups and agencies as to what the objectives of u~ban planning, both in general and in particular, are.

(c)

The paucity and complexity of the models of system behaviour

on ~hich

good optimisation models must be founded. (d)

The need to treat the system in a hi@h disaggregated form to get anywhere near to reality - this produces data and computational problems of a high order.

(2)

0ptimisation techniques have been relatively little used in urban planning in general, although important applications have been made in some narrowly defined areas.

This is largely because of a lack of sympathy between system

scientists and urban planners brought about by a failure by systems sciautists on the one hand to set their work in a context seen to be relevant to urban planners and on the other an unreasonable mistrust by many urban planners of quantitative techniques.

(3)

We are now at a stage where it is possible to construct general models of urban systems which have a reasonable degree of correspondence with the concepts of urban planners and one sufficiently formally defined to allow meaningful optimisation within and in some instances between individual sectors.

However

the value of sectoral optimisation models will be enhanced by being set in an overall analytical framework.

Beware excessively behavioural assumptions as

these have prevented much of the interesting work in regional science and urban econometrics from finding its way into urban planning practice.

(4)

Finally both the challenge and opportunities (see Figures lO and ll) are greater than ever before.

The growing problems of cities are ~here for all to

see and yet this is one of the aspects of modern society where modern problem solving methods have made least impact.

Even if systems analysis and

optimisation techniques do no more than provide a formal arbitrator in the plan development process they will have done urban planning a great service. FIGURE l0 CH3T,!,~GE (1)

Cities the cradle, haven and tomb? of civilisation.

(2)

Urban problems of an unprecedented kind and degree.

(3) Modern problem solving techniques little applied. FIGURE ll OPPORTUNITIES (1)

Increased availability of the right kind of info~nation.

(2)

Better understanding and behavioural models.

(5)

Continuous improvement in problem solving techniques.

(4)

Developing sympathy amongst urban planners.

SELECTION OF OPTIMAL INDUSTRIAL CLUSTERS FOR REGIONAL DEVELOPMENT Stan Czamanski I

The Research Problem The importance of a new industry to a depressed region resides not only in the volume of new employment and income which it provides but very often primarily in its indirect or nmltiplier effects.

A common feature of depressed regions is the

general weakness of multiplier effects capable of being generated in their economies. It is due mainly to the patchiness of interindustry linkages.

The absence of sub-

stantial indirect effects which ordinarily accompany new investments constitutes one of the greatest obstacles to efforts aimed at invigorating the economies of underdeveloped regions.

The weakness of the multiplier effects makes remedial programs

expensive and the resulting progress slow. 2 The problem can also be viewed in a slightly different way.

The study of indus-

trial location patterns brings to the fore the phenomenon of significant progressive clustering of economic activities in a small number of urban-industrial agglomerations.

This remarkable feature of locational preferences of industries is increas-

ingly exploited for fostering regional progress by promoting the emergence of growth poles or more generally by furthering spatially imbalanced development.

Explanations

of the emergence of spatial concentrations of industries revolve around the extent to which geographical proximity between certain classes of plants confers significant advantages.

The benefits attending the reduction in the friction of space may be due

to savings on transportation costs, especially in the case of products which are weight-losing, transported in hot state~ or capable of being transferred over a short distance without packaging, by pipes, belts, or conveyors.

In other cases, especi-

ally those in which storage and related interest costs are substantial, and advanced planning difficult, the advantages associated with spatial proximity may be due to savings in transportation time rather than in transportation costs.

Many industries

are attracted to existing clusters because of the importance of human face-to-face contacts~ of presence of external services, of an existing pool of trained labor, or more generally because of the possibility of realizing savings either in investments or production costs. In a depressed, open regional economy a major breakthrough can result only from iThe author is Professor of City and Regional Planning at Cornell University, Ithaca, N. Y., U. S. A., and Research Associate at the Institute of Public Affairs, Dalhousie University, Halifax, N. S., Canada. The financial support of the Economic Development Administration, U. S. Department of Commerce, and of the National Science Foundation is gratefully acknowledged, as is the research assistance of Abby J. Cohen and Stephen B. Ellis. 2For some early studies bearing on this problem see [6, 5, 4]; and for later work.[~].

10 the introduction of new productive activities.

Hence the empirical question fre-

quently faced by regional planners or political decision makers is the selection of industries to be supported or attracted into the region.

Recognition of the impor-

tance of the multiplier effects for both short-run and long-run development strategies adds a new criterion to those customarily used, such as output-labor, capitallabor, or capital-output ratios, relative rates of growth of industries considered, or the extent of required outlays on infl-astructu~e. Increased interest in the theor5, of imbalanced ~ o w t h and in the ~ o w t h poles hypothesis has led to the fommulation more recently of the idea of promoting clusters of related industries rather than scattered individual plants. 3

The main rea-

son seems to be a desime to realize some of the external economies which are internalized in a cluster of related activities.

Moreover, it is assumed, so far without

any empirical proof, that the ~ o w t h of multiplier effects in a region experiencing an influx of new activities is more than proportional to the gmowth of its economy. The hypothesis is based on the plausible idea that the introduction of new indus~ t-Dies pro~essively reduces leakages and reinforces the indirect impact of new activities until a point is reached when in a group of related industries linked by flows of goods and services the multiplier effects become significantly stronger, signalling a qualitative breakthrough. It is, however, not immediately obvious which industries form groupings subject to external economies.

Traditionally, industrial complexes have been identified on

the basis of studies of either technological processes, or of spatial associations. Both methods present some distinct advantages but also some rather obvious difficulties. 4 A different approach was followed in the reseamch reported here, with major emphasis focussed on subsystems of industries capable of generating external economies and identified on the basis of flows of goods and services connecting them either directly or often only ~ndirect!y. 2.

Identification of Clusters of Industries The main data basis for the study was the United States 478 x 478 input-output

table for 1963 reduced to size 172 x 172. 5

The grouping of the matrix was carried 6

out with the help of a program described elsewhere.

3For an overview see [13]; and for some of the original points of view, [16] and [2].

4The various approaches ame presented in [14; 8; I0; ii; 7; 15; 17; and 1]. 5The reduction was necessary because such detailed data could not be reconciled with those available for the study of spatial patterns for which SMSA's were the basic units. In addition the operation of such a large table presented numerous computing difficulties involving very high costs for computem time which did not appear to be justifiable. 6For details see [12].

Clusters of industries were defined as groups of sectors with relatively stronger ties among themselves than with the rest of the economy.

Three different but

complementary methods were developed for the purpose of identifying clusters according to this definition.

The first two methods were based on network analysis, while

the third applied principal components analysis.

The first method used as the cri-

terion for including an industry into a cluster a single strongest link between any of the industries already in the cluster and any of the remaining industries.

The

criterion of the second method was the strength of the links with all the industries already in the cluster.

The third method used as criterion for including an indus-

try in a cluster similarity between its total profile of suppliers and customers (not o n ~

those belonging to a cluster).

The research started with an examination

of the relative importance of flows between pairs of industries.

The following four

coefficients derived from the input-output flow table described the relative importance of the links, either for the supplying or for the receiving sector: X..

aij

=

13 Z x.. i i]

bij

=

~] I xij j

=

;

aji

;

b.. 31

X.. ]i

Z x.. j ]i

X..

X..

=

J~ I x.. i ]i

;

where = yearly flow between industry i and j. ~] The starting point of the first method was a triangular E matrix, the elements

X,.

e.. of which were formed by defining aS e.. l] = max (ai~'3 aJ i' bij ' bJ i) for i > j, and e..

=

0 for i ~ j.

The model made use of an e. column vector (i = 1 .... 172) for any of the 172 1 industries (all were used in turn in the starting position). The entries in this column vector were ranked by interchanging rows and columns in the triangular matrix E.

The second entry in the column vector was then the strongest link that the orig-

inal industry (the first entry) had with any other sector in the economy.

This

second industry was represented by the adjoining column vector. Next, the entries in the second vector were ordered and the top entries in both vectors were examined for the strongest link with any other sector in the economy.

In the first variant the process was repeated until the mean value of the

entries among the industries in the cluster, or ~ ~ eij/n for i, j = l...n, began to decline, where n is the number of industries in the

cluster.

The second method also makes use of the triangular E matrix defined above, and follows a similar procedume of ranking industries by interchanging rows and columns. The ranking is done, however, on the basis of the sum, for all industries already in the cluster, of the links to the other sectors.

This method employed the same stop-

~2

ping criterion that was used in the first method. Not unexpectedly, the clusters identified by the second method were larger than those identified by the first method.

The second method included industries which

have links to several industries in the cluster but no str~ng links to any of the industries in the cluster. Both methods had a tendency to become "side-tracked" whenever industries being added to the new cluster had stronger ties with the industries of another cluster. In such cases the program kept adding to the cluster being formed all members of the tight cluster already identified.

The "side-tracking" would limit the number of

clusters identified to a small number of the strongest groupings.

This had been

expected since the program, while not bent upon finding the maximum maximorum, cleamly tries to find the strongest groupings by including industries with links represented by the highest e.. values. Once industries belonging to a previously identil] fled tight cluster were included in the cluster being formed, and the mean value of links began to decline, the program stopped without revealing the weaker grouping of possibly great interest. In order to overcome this difficulty, the preliminary clusters defined by the fimst and second methods were supplemented by including all industries with links (e.. coefficients) greater than .200. This additional procedure was carried out 13 with the help of a computer program which created for each of the 172 industries lists of all significant e.. l] coefficients in decreasing order. Hence, to be included in the cluster an industry had to meet at least one of the following criteria: (a)

More than two-tenths (.200) of its output was absorbed by one of the industries of the preliminary cluster,

(b)

More than two-tenths (.200) of the output of one of the industries in the preliminary cluster was absorbed by the industry examined,

(c)

More than two-tenths (,200) of inputs of the industry examined came from one of the industries in the preliminary cluster, or

(d)

More than two-tenths (.200) of the inputs of one of the industries in the preliminamy cluster came from the industry examined.

The analytically more interesting indirect links among industries were identified with the help of the third method~ described elsewhere. 7

It also started with

the 172 x 172 input-output flow matrix, from which an n x #n matrix of zero order correlation coefficients was derived. r

=

[r(aik.ail)

I

r(bki.bli):,

!

r(aik.bli) !

,I!

r(bki.ail)~

;

Next an n x n covariance matrix was for~ed in terms of derivations from the mean values.

K

:

E

(r-

7For details see [9].

)TJ

13

In order to identify, from the set of all industries, the subgroup belonging to a cluster, an iterative process was applied, eliminating all industries having a null column or a null row vector.

The relative strength of the links binding the

remaining industries together was assessed with the help of eigenvalues of the R matrix.

The ratios of the characteristic roots to the trace of the R matrix define

an Index of Association

C n

-

n tr R

x i00

where k

=

eigenvalue, or characteristic root.

This provided an aggregate measure of the strength of the ties connecting the industries remaining in the R matrix - each C. indicating the existence of an industrial cluster.

The indirect links revealed by the third method were analytically

highly significant since two industries k and 1 may be members of an industrial complex in the absence of direct flows between them.

The three methods, using different

criteria, yielded results fully consistent with one another when applied to the Washington State input-output matrix and to the United States 1963 matrix.

The second

method yielded larger but less closely linked groupings of industries than the first, while the application of the third method resulted in still larger and more diffuse clusters. 3.

All three were, however, consistent in their ranking of industries.

Characteristics of Industrial Groupings Seventeen clusters were identified in the 1963 U. S. economy.

Two groups of

industries which were too small and too weakly linked to be defined as clusters are also of some interest.

Characteristics of the seventeen clusters are listed in the

summary tables. The first and most obvious classification of the seventeen clusters is by manufacturing or service, according to the nature of industries of Which they are composed.

There are few service sectors in the predominently manufacturing clusters,

and vice versa.

The general weakness, or rather relative unimportance, of links

between manufacturing and service sectors was unexpected.

It runs counter to any

notions of complementarity and may have significant implications for regional development theories and strategies. The seventeen clusters differ greatly, in terms of the number of industries included, and in terms of the total size of the cluster as measured by value added. The two smallest clusters, Recreation and Government, are each composed of only eight sectors, while the largest cluster, Construction, contains forty-two sectors. More striking is the range in terms of value added for each cluster.

The two largest

clusters, in terms of output, are Real Estate and Construction, with 171.6 and 155.9 billions of dollars of value added respectively.

The smallest is the Services

cluster, which has a value added figure of 25.0 billion dollars, Table i. Examination of the identified clusters revealed that the component industries

14

widely differed in their relative importance to the cluster as a whole.

While the

removal of some of the sectors would have little effect upon the structare of the cluster, in many clusters the exclusion of a single industry would lead to disintegration.

These important sectors were called central industries and were defined as

having at least four links to other industries in the cluster.

The centrality of an

industry may be measured not only by its degree (number of links), but also by the strength of the links. TABLE 1 SIZE OF CLUSTERS

Cluster 1 2 3 4 5 6 7 8 9 i0 ii 12 13 14 15 16 17

Number of Industries

Size in $000 of Value Added

21 42 12 15 15 22 l0 13 23 10 14 18 23 12 8 ii 8

54,823,590 155,886,704 40,527,318 45,349,988 34,321,747 71,411,137 38,632,704 56,753,656 60,501,678 28,503,745 47,534,129 55,240,961 171,612,897 24,973,159 50,034,347 81,237,289

Foods & Agricultural Products Construction Textiles Wood & Wood Products Paper $ Printing Petrochemicals Petroleum Leather Products Iron & Steel Nonferrous metals Communications g Electronics Automotive Real Estate Services Recreation Medical Services Government

34,136,089

Percentage of GNP* 9.39 26.59 6.94 7.76 5.87 12.23 6.61 9.71 10.36 4.88 8.14

9.46 29.39 4.27 8.58 13.91 5.84

~Percentages of GNP are not additive because many sectors are members of more than one cluster. Each of the seventeen clusters contained at least one central industry; some clusters had two or three.

It may be noted that nearly all of the central indus-

tries appear to possess some special technological significance in their respective clusters.

This is a reflection of their position in the sequence of processing op-

erations.

Interestingly, however, central industries are not always the largest ones

in a cluster.

In only ten of the seventeen clusters were the central industries the

leaders in terms of output or value added, and even then they were not significantly larger than some other industries. The importance of a central industry can be measured in a number of ways, as shown in Table 2.

The measures used were size of the industry, the number of links

to other industries (degree), the percentage of all links to the central industry, and the relative average strength of links. A simple way of classifying the clusters would be to distinguish between single and multi-centered clusters.

Of the seventeen clusters considered eight groups had

Foods ~ Agricultural Products Construction

Textiles

Wood & Wood Products Paper & Printing

Petrochemicals

Petroleum Leather Products

Iron & Steel Nonferrous Metals

Communications & Electronics Automotive Real Estate Services

Recreation

Medical Services

Government

1

2

3

4

5

6

7 8

9 l0

ii

15

16

17

12 13 14

(2)

Name of Cluster

(I)

Cluster N~mber

7,708 1,819 2,053

4 4 5

5

8

13,710 2,904

12 5 34 4 4 4 8 6 8 4 6 7 8 4 4 6 4 4 16 6 5 4 5 15 17

(5)

22,104 1,683 34,830 786 395 6,679 3,540 520 1,573 2,961 1,857 6,171 2,865 I~322 3~358 273 91 33 8,617 1,013 2,128 5,341 7,867 12,781 29,759

(4)

Size of C.I. in $i,000 Value Added Degree

Agriculture & Related Services (i) Nonmetallic Mineral Mining (8) Contract Construction (i0) Hydraulic Cement (84) Floor Covering Mills (13) Apparel (14) Fabric g Yarm Mills (30) Logging (35) Sawmills (36) Commercial Printing (46) Paper Mills (48) Basic Chemicals (60) Fibers, Plastics, Rubbers (61) Tires & Inner Tubes (70) Petroleum Refining (68) Leather Tanning, Finishing (75) Boot, Shoe Cut Stock (77) Leather Goods, n.e.c. (82) Steel Rolling, Finishing (90) Primary Nonferrous Metals (92) Nonferrous Rolling & Drawing (94) Communication Equipment (120) Aircraft ~ Parts (124) Motor Vehicles & Equipment (123) Real Estate (158) Miscellaneous Business, Personal & Repair Services (160) Amusement & Recreation Services, n.e.c. (166) Physicians, Surgeons, Dentists, etc. (167) Medical & Health Services Government Enterprises (171)

(3)

Central Industries

TABLE 2 CHARACTERISTICS OF CLUSTERS

28.8 28.8 62.5

71.4

72.7

54.5 10.6 72.3 8.5 23.5 23.5 47.0 35.3 35.3 25.0 37.5 29.2 33.3 16.7 40.0 37.5 25.0 25.0 72.7 54.5 45.5 25.0 31.3 88.2 68.0

(6)

Percentage of Links to Central Industry

.448 .443 .461

.452

.397

.402 .461 .461 .461 .383 .383 .383 .386 .386 .431 .431 .369 .369 .369 .560 .454 .454 .454 .417 .371 .371 .391 .391 .280 .353

(7)

Average Strength of All Links

.429 .380 .569

.387

.425

.451 .288 .524 .373 .368 .559 .358 .370 .385 .453 .461 .345 .348 .398 .561 .455 .545 .343 .405 .475 .312 .427 .474 .280 .344

(8)

Average Strength of Links to Central Industry

0.968 0.858 1.234

0.856

1.070

1.121 0.625 1.136 0.809 0.961 1.460 0.934 0.959 0.997 1.051 1.070 0.934 0.943 1.065 1.001 1.002 1.200 0.756 0.971 1.280 0.841 1.092 1.212 1.000 0.974

(9)

Ratio 8 ÷ 7

16

one centTal industry, five had two, and four clusters had three central industries. However, the number of central industries is not neamly as important as the relative strength of these industTies.

Some of the clusters appeared To be dominated by a

single industmy, while the influence of central industries in other clusters was not nea~ly as great.

The degree of dominance can be judged by the number and intensity

of the links involving the central industry.

The important indicators of dominance

appear in columns 5, 6, and 9 of Table 2. Another measure of the tightness of a cluster is the mean strength of interindustry links, or ~ ~ ei~/n, where n is the number of industries in the cluster. J

Considering only those links for which e.. k .200, the mean value of links ranged l] from .280 for the Automotive cluster to .560 for the Petroleum cluster, indicating that, on average, the industries of the latter were the more interdependent, in terms of interindustmy flows. In order to proceed with the analysis of the structure of the various clusters, a diagram was constructed of each identified grouping.

The diagrams are two-

dimensional representations of the multi-dimensional clusters.

Each circle repre-

sents one sector; central industries are identified by a double circle.

The area of

each circle is proportional to the size of the industry it represents, using value of shipments as an index of industry size.

Value added has been used in the cases

where value of shipments cannot be obtained, or is meaningless, as for some service industries. A line between two circles represents a link between two industries.

The length

of a line between a given pair of industries is inversely proportional to the respective e.. coefficient. In cases where an industry is a member of another cluster, 13 the industmyappeams within a dotted box indicating that there are a number of links to the other cluster that are not shown. Among the seventeen clusters there were five yielding tree diagrams.

In these

clusters there is only one link, om set of links, connecting two industries.

The

clusters giving tree diagrams were Iron & Steel, Automotive, Petroleum, Recreation, and Services.

Clusters which did not form trees, containing two or more circuits

between industries, were called complex groupings.

A small subgroup of clusters has

been termed "snowflakes" because of the appearance of their diagrams.

To be classi-

fied as a "snowflake~ a cluster had to meet the following three conditions:

(i)

have only one central industry, (2) be classified as a tree, and (S) have at least two-thir~s of all links in the cluster involve the central industry. clusters were identified. Services.

Four such

These were Iron G Steel, Automotive, Recreation, and

Iron & Steel proved to be a typical representative of this group, of

great interest for a number of reasons.

The sequence of operations based upon engin-

eering considerations is well preserved

or only slightly distorted by the use of

SIC classification. ing.

It is dominated by its central industry, Steel Rolling & Finish-

The strong revealed ties to the Construction and Automotive clusters are not

unexpected.

17

U.S. E C O N O M Y

II:~ON AND STEEL CLUSTER 3

Iron & ferroalloy ores

6

Authracite, lignite & bituminous CO~I mining

II

Oil & gas field services

23

Canned & frozen foods

27

Beverage industries

h2

Office fu~ibure

44

Partitions & fixtures

90

Steel rolling & finishing

91

Iron & steel f o ~ r i e s

96

Frimary metal industries, n.e.c.

97

metal cans

98

Cutlery, hand tools, hardware

iO0

Fabricated structural metal products

lO1

Screw machine prods. & bolt

102

Metal stamp~ ags

103

Coating, plating, polishing engraving

104

Fabricated wire products, n.e.c.

105

Fabricated metal products, n.e,c.

125

.....................

Ships & boats*

126

Railroad equipment

145

Water transportation*

150

Electric, gas & sanitary service*

158

Real Estate*

___$ *Industry size given by Value Added rather than Value of Shipments.

The Automotive cluster is itself highly significant, not so much in terms of its structure as because of its economic importance, its strong links to other groupings, and the presence of several nodal industries. Of considerable significance for regional development strategies may prove to be a small gro~@ of industries, called nodal industries, which belong to more than one cluster or have strong links to another cluster.

One industry, Contract Con-

struction, was included in eleven of the seventeen clusters. and Automotive~ each belonged to six clusters.

Two others, Real Estate

The following is an enumeration of

the nodal industries and their number of occurrences: Industry Contract Construction Real Estate MotoP Vehicles Agricultural Products Advertising Services FlooP Covering Mills 13 Others 58 Others

No. of Clusters Appeared In ii

6 6 5 4 4 4 3

2

All of the nodal industries listed above are also central industries in their respective clusters, and most were highly correlated with population.

The presence

18

of these nodal industries is apparently highly significant for the formation of spatial industrial complexes. KK

ECONOMY

r ........

AUTOMOTIVE

CLUSTER

f

I ......... NICATIONS

~ ~ T R ~ N I C S

kCONSTRUC]IC~

I

"r I

'I J

I

I

!

'

Io'~

]

IcLus+E~ ITEXBLESI

96

...... "~

~

Primary m e t a l

industries,

Cutlery,hand tools,hardwere

101 I02 109 114

Screwmachineprods. & bolts Metalstamplnge Metalworkingmachinery Machinery~exceptelectrical,

119 122 123 124 130 16~

Radio,TV receivingequip. Electricalprods.,n.e.e. M~torvehicles& equip. Aircraft& parts Mechanlealmeasttrin~ devices Automobilerepair & service*

n.e.c.

l --ti

L.. . . .

Contractconstruction Floorcoveringmills Fabricatedtextiles,n.e.¢. Tires& i n n e r tubes Iron& steelfoundries Nonferrousfoundries n.e.e.

~2 ---------'O

v

I0 13 34 70 91 95

------i

I ! i

I

IPE+~C E~+ t

~NONFERROUS I

i

IM£TALS

hRoN&

S'~EEL CLUSIER I- . . . . . . . . . . . . . . . . . . . . .

* Industrysize givenby Value Added ratherthan Valueof Shipment8.

I

L~us_L~_ J

Several Of the identified clusters are worthy of special consideration, either because of their structure or their relevance in any planning strategy. most important groupings was the Construction cluster. its ubiquitous nature,

One of the

This was due to its size and

The central industry of the cluster, Contract Construction,

repeatedly appeared in other clusters and made a substantial contribution to the value added totals of these clusters,

Diagramatically, the Construction cluster was

very nearly a pemfect "snowflake", except for the presence of a small subgroup centered around the Nonmetallic Mineral Mining industry,

The weakness of the ties

within this subgroup did not justify defining it as a cluster, and it remained a part of the Construction group. The Real Estate cluster was similar to Construction in two important respects. First, its central industry was extremely large and appeared in many clusters.

Sec-

ondly, the diagram of the cluster had a "snowflake" appearance, but the presence of circuits prevented its classification as such. Several clusters, particularly those primarily engaged in manufacturing rather than providing services, were composed almost exclusively of technically related industries.

An example of this is the Textiles cluster, although the tree-like se-

quence of processes is obscured by the circuits present in the diagram.

The lack of

19

CONSTRUCTION CLUSTER

~S, E C O N O M Y

IME~ALS

6

CLUSTER !

8 i0 13 19 36 37 39 41 43 44 51 6h 69 79 8~ 85 86 87 88 89 90 94 99 lOO

TEXTILES

U.S. E C O N O M Y

CLUSTER

104 108 113 115 117 118 119 123 124 130 142 147 150 158 160 163

ICLUS~£R 171

I I

C i LOS~ER~' ~ I

!

:

Anthracite, li@nite & bituminous coal mining Nor~metallic m~neral mining Contract construction Floor covering mills Scroll arms ammunition Sawmills & planing mills Millworks & relate~ product~ Misc. wood prods. & wood finishing Household furniture Public bldg. & related furniture Partitions & fixtures Bldg. paper & board mills Paints, varnishes & allied products. Paving & roofing materials. Leather gloves & mittens Hydraulic cement Structural clay products Pottex~ & related products Concrete, ~/ps%~ & plaster products Cut stone & stone products Nor~netallic mineral prods., n.e.c. Steel rolling & finishing Nonferrous rolling ~& drawing Plumbing & nonelectric heatiag equip. Fabricated structural metal products Fabricated wire prods., n.e.c. Construction & llke equ~,p. Service industry machines Electric distribution prods. Household appliances Lighting & wiring devices Radio, TV receiving equip. Motor vehicles & equip. Aircraft & parts Mechanical measuring devices Railroads & related services* Pipe line transportation* Electric, gas & sanitary ser.* Real estate* Services: misc. business, mist repair, personal* Prof. services (except education), non-profit organizatians* Government enterprises*

*Industry size given by Value Added rather than Value of Shipments,

i

iO

Contract construction

13

Floor covering mills

14

Apparel (except~

30

Fabric & yarn mills ; textile flulshln 6

31

Narrow fabric mills

33

Textile gOOdS~ n.e.e.

goeds)

Knitting mills

i IAU ~OMOYIVE

34

Fabricated textiles, n.e.c.

61

Fibers, plastics, rubbers

71

Rubber footwear

123

Motor vehicles & equip.

I~O

Costume Jewelry & notions

20

an orderly transition from one operation to the next is partly due to SIC classifications and the grouping required to obtain spatial data.

In one of the central in-

dustries in this cluster, for example, the operations of spinning, weaving, and finishing have been lumped together in one sector.

Moreover s the cluster is comprised

of industries using different raw materials to satisfy heterogeneous needs, both natural fibers and synthetics being included in the cluster. One surprising finding was the absence of significant links between the Petroleum and the Petrochemicals clusters.

The links between them are technically cru-

cial but not economically significant.

The Petroleum cluster has an almost tree-

like structure, with Petroleum Refining occupying the central position.

The se-

quence of the various operations is reflected in the diagram of the cluster.

The

Petrochemicals cluster only weakly exhibits the expected sequence of technical process.

The tree-like structure familiar from descriptions of petrochemical complexes

is lost, at least partly because of the broad groupings of operations found in the SIC classifications.

The cluster is unusual in that it could be easily cleaved into

two separate groups, if not for the extremely strong link between the Basic Chemicals industry and the Fibers, Plastics, and Rubbers industry. U.S. E C O N O M Y

RETROCHEMICALS CLUSTER

~!

LEc~ND Agricultural services, huuting, trapping; fizharies* Chemical & fertilizer mineral mining

72

lO ( ~

A3

Contract construction Floor covering mills

30

Fabric, yarn mills; textile finishing

31

Narrow fabric m i l l s

33 6O

Textile goods, n.e.c,

61

Fibers, plastics, rubbers Paints, varnishes, and allied 1 ~ u c t s

64 65

Basic chemicals

GI~ & wood chemicals

Agricultural chemicals Misc. chemleal products

[ AGRICULI~RI 70

T i r e s & inner tubes Rubber footwear

Reclaimed rubber

9

Fabricated rubber producte~ n.e.e. Misc. plastic products S~eci~l industry machinery

V

Mortar vehicles & equip.

67,

i I .L

131~ 141

Photographic e ~ p . & supplies Misc. manufacturers

*Industry size given by Value Added rather than Value of Shipments.

The above analysis provides insights into the extent and nature of intersectoral interdependence.

However, it is unclear whether the flow of goods and services

21

between s e c t o r s can help explain the phenomenon of spatial aggTegation.

Obviously,

the relative value of flows between two industries need not be related to their mutual attraction in space.

Work now under way is attempting to determine the cor-

relation between flows and spatial association. REFERENCES

[l] Bergsman=, J., P. Greenston, and R. Healy~ "The Agglomeration Process in Urban Growth", Urban Studies, IX, 3 (1972).

[2]

Boudeville, J. R.

[3]

Casetti, Emilio, "Optimal Interregional Investment Transfers", Journal of Regional Science, Vlll, 1 (1968).

[4]

Chenery, Hollis B., "Development Policies for Southern Italy", Quarterly Journal ofEconpmics, LXXVI (November 1962).

[5]

Chenery, Hollis B. and Paul G. Clark.

[6]

Chenery, Hollis B., Paul G. Clamk, and V. Cao-Pinna. o f the ' ItalianEconomy, 1953.

[7]

Collida, A., P. L. Fano, and M. D'Ambrosio (F. Angeli, ed.). mico e Crescita Urbana en Italia, 1968.

Problems of Regional Economic Planning, 1966.

Inter industry Economics , 1962. The Structure and Growth Sviluppo Econo-

[8] Czamanski, Stan, "Industrial Location and Urban Growth", The Town Planning Revi__~e__w,Liverpool, XXXVI, 3 (October 1965).

[9]

• "Linkages Between Industries in Urban-Regional Complexes", in G. G. Judge and T. Takayama, Studies in Economic P!anning Over Space and Time, 19"73.

[1o]

.

[ll]

_ .

~'A Method of Forecasting Metropolitan Growth by Means of Distributed Lags Analysis"~ Journal of Regional Science, VI (1965). "A Model of Urban Growth", Papers, Regional Science Association, XlIl (1965).

[12]

Czamanski, Stan, with the assistance of Emil E. Malizia, "Applzcability" and Limitations in the Use of National Input-Output Tables for Regional Studies", Papers , Regional , Science Association, XXIiI (1969).

[13]

Hez~ansen, T. "Development Poles and Related Theories", in N. M. Hanson (ed.) G_r,owth Centers in Regiona ! Economic Development, 1972.

[14]

Isard, Walter, Eugene W. Schooler, and Thomas Vietorisz. Analysis and Regional Development, 1959.

[15]

K!aassen, L. H.

[16]

Perroux, F. "Economic Space: Economics , LXIV (1950).

[17]

van Wickeren, A. "An Attraction Analysis for the Asturian Economy", Regional and Urban Economics, II, 3 (1972).

Industrial Complex

Methods of Selecting Industries for Depressed Areas, 1967. Theory and Applications", Quarterly Journal of

OPTIMAL INVESTMENT POLICIES IN TRANSPORTATION NETWORKS

(~) S. Giulianelli - A. La Bella

i. INTRODUCTION A relevant feature of land use planning is the design of transportation networks. The problem which public administrators are often confronted with is that of apportioning a limited budget to the various branches of a network so that to achieve the goal put forward by the administration. Namely, given the demand for transportation among the centers connected by the network, one may wish to find the optimal investment policy, in order to minimize the total transportation time. The main difficulty of this kind of problems is that the objective function turns out to be neither convex nor separable whereas the constraints turn out to be non linear. Some authors have tackled the problem using, in our opinion, too many sim plifying assumptions. For example, the transit time on each branch of the network is assumed to be independent of the flow (I) (2) (3). Moreover, in case of several origins and destinations only zero-one investments in each branch are considered (I) (2), whereas discrete investments are analyzed for one origin only (3). Although the hypothesis that the transit time is independent of the flow seems to be quite coarse, we keep it in this paper, where, as our first contribution to the problem, we propose to overcome all other restrictions. Namely, we consider the case of several origins and destinations with the investments in each branch taken as continuous variables. We also suggest a heuristic computational procedure which appears to he more efficient compared to those so far appeared in the literature.

2. DEFINITION OF THE PROBLEM The basic hypotheses of this work can be summarized as follows:

Hypothesis 1: (transit time hypothesis): the transit time tZ on each arc Z is indepe~ dent of the flow.

Hypothesis 2 (behavioural hypothesis): the traffic flow distribution is such that the overall transportation time is minimized.

Hypothesis 3: the dependence of the transit time t% associated with arc % on the apportionment cz is of the following linear type with saturation

(~) Sandro Giulianelli - Agostino La Belle - Centro di Studio dei Sistemi di ControlIo e Calcolo Automatici, C.N.R. - Istituto di Automatica, Unlversit~ degli Studi di R o m a - Via Eudossiena, 18 - 00184 Rome.

23

where ~,

~£, c ~

are non negative constants.

Then the problem can be stated as follows:

Given: 1 - A transportation network G(N,A), where N is the set of nodes and A the set arcs.

of

2 - The transportation demand matrix R = {r..}

13

3 - A budget B 4 - An objective function T

T =

Z

f(~)~j~ij

(c)

(I)

ieN jeN (i.e. the overall transit time), where k k x.. : column vector, whose components xij represent the flow on the k-th path Pij --zj between i and j, i,jeN due to the demand rij ; k=l,2,...,qij qij : number of different paths between the nodes i and j, i,jeN P~. : set of branches in the path p~j A_ij : incidence matrix arcs-paths for the origin- destination pair i-j, i,jeN, whose entries am,k, meA, k=l,2 .... ,qiJ , are defined by I I

for

0

for

am'k =

m e e~. 13 m ~ P,. k ij

(2)

~(c): column vector with components tm(Cm) , meA : column vector with components Cm, meA Minimize T, subject to the following constraints qij

k

X

~,,(c) k= I I~ - -

Z

< B

c m

= r.. 1J

VieN, YjeN

(3)

(4)

-

mcA c

> o

(5)

e

< c

(6)

24

where C -- {c } , --

meA

m

3. PROPERTIES OF THE OPTIMAL SOLUTION For each arc meA, it is possible to define a return as =mYm, where Ym is the total flow on the arc. Then one may think to obtain the optimal solution by putting the branches in o[ der of decreasing return and investing in this order. Actually, this approach would be fallacious, because by investing on each branch we change the relative transit time. Consequently a new distribution of flow may follow: this phenomenon modifies the problem due to the possible variation of the return associated w~th each branch. It is evident that the investment on some branches albeit optimal relative to a certain amount of capital may fail to indicate the optimal solution for a larger bu~ get. Nevertheless the following properties and theorems hold: PPop~rty 1 - The whole traffic demand for each origin- destination pair runs along one and only one path. In fact from the hypotheses I) and 2) it results that each individual chooses the minimal time path from the origin to the destination (i.e. there is no interaction among individuals). R e ~ P k 1 - For each investment vector c the problem of determining the traffic distribution in the network is reduced to that of finding the shortest route from each origin to all destinations. This immediately follows from Property I. PPoperty 2 - The flow originating at each origin branches out into a tree. Therefore, for each vector c the total traffic flow is distributed on a network determined by the union of as many trees as these are origins. If now we put:

X(£) = I

[ ~j ~ij(£)

ioN jeN where ~(~) is a column vector with components Ym(C), msA, the following theorem hold: THEOREM i - Suppose that the optimal flow distribution [e(B) is known. Consider the subnetwork H e = G(NH~,AHe) obtained from G(N,A) cancelling out all branches where the flow is zero (cfr. Property 2). Then the optimal investment policy c~(B) is obtained by investing the maximum amount of capital c"~ in each branch in order of decreasing ~eturn until the budget is reached. Proof - Under the hypotheses set forth, the problem (I), (3)-(6), is reduced to the following linear programming problem: min

[ £eA H

W

y£(B£

e£ c£ )

25

< B EAH, whose optimal solution is the investment policy stated in the theorem.

4. EXHAUSTIVE PROCEDURE Based on theprevious sections, the following procedure leads to the optimal investment policy. Step 1 - Construct the set H of subnetworks H = G(NH,A ) combining in all possible ways one tree for each origin to all destinations for ~hich there is demand for transportation. Assume that the traffic demand is satisfied using the tree chosen for the relative origin only. Let [(H) be the vector of traffic distribution. Step

2 -

min

For each HeH

[

solve the following linear prograu~ing problem:

Ym(H)tm(Cm ) =

m~AH

I

Ym(H)tm[C~ (H)]

mcAH

c < B m meA H (for the solution procedure cfr. theorem i). Step

~ -

Compute

min [ Ym(H) tm Item(H)] HeH msA H Con~nent - The exhaustive procedure singles out a finite number of points from the subset of 1~he Euclidean space whose coordinates are Cm, meA, defined by

I

< B

c m

-

m~A

c

> 0 m

5.

¥ mcA

-

HEURISTIC PROCEDURE

The procedure previously described examines a finite and discrete subset of solutions belonging to the infinite and continuous set of possible solutions, and chooses among them the optimal one. Now this subset of admissible solutions, already defined, can be made up of a very large number of elements. It is therefore desirable to obtain an algorithm that makes it possible to find the optimal solution more rapidly. This algorithm can be set forth as follows:

26

Step I - Find the optimal traffic distribution, corresponding to the initial s! tuation of no investment. Invest on the arcs in order of decreasing return until a new distribution of traffic is obtained. Step i - Cancel the preceeding investment and, assuming traffic distribution ob, tained in the previous step, invest on the arcs in order of decreasing return u~ til a distribution of traffic is obtained for an amount of capital greater than the previous investment. If no redistribution is obtained, then stop. The optimal solution is the one that corresponds to the investment made at the K-th iteration, such that: Tk(B) = min T.(B)I i where Ti(B) is the overall transit time obtained in the i-th iteration for an investment equal to B. Notice that some traffic distributions may occur in the euristic procedure which are not stable for the value B of the budget. In such cases, in order to comp~ te Tk(B), these distributions are replaced by those obtained starting with the unst~ ble distributions with an investment equal to B.

6. NUMERICAL EXAMPLE Let us consider the network of fig. I, with the following traffic demand matrix:

L131

R=

0

0

4

0

0

5

3

6

Fig. I Let the dependence of the transit time on each arc on the investment be as follows: ta = I 0 - 2 c t

a

b = 6-c b

C

= 6-c C

O

to one

The different

here is the following.

minimal

Z

combining

(1971).

2. MEDIUM OPTIMAL

i.e.

either

out how much has to be done even about the best known

problems

~

may

algorithm.

help to point

There

algorithm - Van Slike

with respect

on the other.

but suitable

with

of the tree will be

per arc,

factor or optimizing

some constraint

Most of this work could be considered

cij

total weight

~ Pij

Assume

we have

112

and o b v i o u s l y the guessed value of z is too pessimistic

since there

exists a tree T for w h i c h a smaller value can be found. w

>

0

Similarly Z Z


n-1

...

S

f

(f K)

(2.17) D

Theorem 2. 3. An n-step solution ~n yields one and only one sequence of states, S [~]

, s.t.:

S[%n] is called n-step solution sequence, and it is composed of n-1 states called intermediate solution states. n P ~ o f . If ~ is am n-step solution, because of Definition 2.5 and relations(2.14) and (2.17), it individuates a sequence of n-1 states; the sequence is unique since in (2.17) each a is, because of (2.15), an operator Yi' and, therefore, because of Definition ~. I, each a is a function.

D We can use the mathematical notion of directed graph in order to handle the SSPS, both in its aspects of representation and search. We associate a vertex to each state and an arc to each operator. -

When we apply an operator "(i to a state s, we obtain a new state s* = In the gr_aph model we shall have an which joins s and s* .

arc (corresponding to the operator yi )

If we define spree "cost___~s" associated with the operators, we could be interested in obtaining an optimum solution", i.e., a solution whose total cost, i.e., the sum of the costs of the a j & ~ , is minimal. More precisely~ we may introduce the following definitions. Definition 2.6. An operator yi~ F is associated with a cos__~tCi, which evaluates the application of Yi between a state s and a state t, i.e.: s

~q--~--*t

( 2 . ~9)

Definition 2.7. The cost of a solution ~ .

%n is defined as :

=

F] (2.2o)

n

J

(2.21) []

150

Definition 2.8. Given a problem P and all its solutions ~ , w e

of the_:solution to P, the nimum of the costs solution ~ A the related n-step solution to P.

call optimum cost

and

call o timum

In the graph model the optimum solution is represented by a minimum cost path from i to f. The search aspects of the automatic problem-solving consist mainly in obtainig a solution, possibly an optimum solution. The construction of the solution (or of the optimum solution) involves the use of some search strategies. Since the dimension of the state space, for real problems, is usually of big magnitude~ the problem arises of limiting the occupation of memory and the compu tation time. Thus, normally, the graph associated to the state space, is not stored in a way that is called explicit description, i.e., with the complete list of all vertices and of all arcs. In effects,the graph is stored in a way,which is denoted implicit description , and which consists of the initial state i,of the set of operators F, and of the final state f (or the set of final states K). Therefore the search strategy always consists of two processes which are developed in parallel : I) incremental construction of the path (i.e., of the solution); 2) incremental explicitation of the graph (i.e., by the application of the n-step global operator F n) . Because of this incremental procedure, these search strategies are called ,ex~ansion ,tehni ques. The most conceptually simple expansion search strategy is the breadth first method~ which consists in computing step by step r n and testing if some of the new states which are obtained is f (or is an alement of K). This strategy is very costly because of the memory and of the computation time which are required. New techniques of expansion are based on the application of F n at each step (i.e., of F2 for n > i) not to all the new states obtained at the previous step, but only to one new state, selected with some criterion (therefore F 2 is practically equivalent to rl). Please note that in this way the expansion process implies an ordering the states which are expanded. We say that a state s is expanded, when we apply r 2 to

s

(or Fi

of

to s).

Depth first method is one expansion search strategy which guides (i.e.which, orders) the expansion process in this way: it is expanded first, the first new state which was obtained in the expansion of the previous step. This method does not assure that the solution which is obtained is an optimum sol~tion. An other expansion technique is the Dijkstra-Dantzig or uniform cost algorithm. In this method a state s is associated with an evalualion function f(s). This function is computed for each state on the basis of the costs of the operators ai which constitute a solution for a problem which has s as final state

151

(and i as initial state). More precisely, f(s) is the cost of the optimum solution for such problem. The value of f(s) for each state s guides the choice of which new state has to be expanded. It is expanded the state s for which is minimum the value

f(s).

Although this method is more efficient than the previous ones~ it is not yet satisfactory for the automatic solution of large scale problems. All these methods use 'blind search strategies" since they are based only on the information which is completely contained inside the representation of the problem. The problem arises of utilizing some additional information, ~alled the "heuristic information", for guiding the search. This information lies outside the given representation and consists of knowledge which is extracted from the "semantic domain" of the problem (Nilsson

(1971)). A classical way of introducing heuristic information, proposed by Hart, Nilsson~ and Raphael, consists in assigning to each state s a new evaluation function ~(s). This function represents an "estimate" of the cost g(s) of a solution to the problem, with f as final state, and with s as an intermediate solution state. More precisely the new evaluation function is the following one; ~(s) = f(s) + ~(s)

(2.22)

where f(s) is the evaluation function of the Dijkstra-Dantzig algorithm and h(s) is a lo_wwerbgund~ of the cost h(s) of the optimum solution for a problem which has s as initial state (and f as final state). Therefore we have :

~(s) _< h(s) ,~(s) < g(s)

(2.23) (2.2.4)

This method, under assumption (2.23), is admissible, i.e., it _nos, whenever it exists, a solution which is an optimum one. An important question is, in this case, how to provide a technique for obtainim4E an evaluation function ~(s) (actually, h(s)). Usually this estimate requires human ingenuity based on an appropriate inspection ~nd processing of the semantic domain of the problem. Therefore the goal of a complete automatization of the problem solving procedure is incompatible with the use of an heuristic information, which is uot contained inside the computer , and which thus implies invention from the man. 9br the above considerations, it seems important to develop new techniques which enable in some way the computer to obtain, within an automatic procedure, the heuristic information itself. This task requires the availability to the computer of new information, out side the SSFS representation~ from which, with an extraction process, the eva-

i52

luation function cau

be computed.

A new way of conceiving the SSPS is therefore required, which we shall call semantic description, in order to distinguish it from the previously exposed syntactic description. The richer information~ which is proper of the semantic description, is embedded on the SSPB representation, by exploding the notion of state, and by associating to each state a structure in which the new inforzation is inserted. We shall describle how to utilize the semantic description in order to obtain the computation of the evaluation function. This technique, which shall be exposed in the following Sections, constitutes a new progress in the direction of automatic problem solving.

III. S~ANTIC DESCRIPTION OF STATE-SPACE APPROACH TO PROBLEM~-SOLVING In this Section we will expose a new formal framework in which it i1~ possible to arrange the new information necessary for the computation of the evaluation function. Since this new formalisation is always related to the SSPS, we will illustra te it with two goals : I. it is sufficient for describing a problem in SSFS 2. it is equivalent to the syntactic description of SSPS. The main idea on which the semantic descriptionis based, is to associate to a state a structure which contains the new information. Therefore we will expose the formal schema which is apted structure of states in SSPS.

to describe the

This formal schema is designed with the purpose of providing these properties: I. it is sufficiently powerful in order to contain a great amount of information from the real problem domain, i.e., the semantics of the problem; 2. it presents great flexibility in order to be utilized for the description of a wide class of problems, arising from very different semantic domains; 3. it is a fruitful basis on which efficient procedures can operate in order to reach the goal of computing the evaluation ftmction, and, in general, the fundamental goals of automatic problem solving. The schema that is proposed here for a formal structure of a state, is quence of attribute and va!ue comples.

a

s___~

For this reason, we shall call the problems, whose structure can be setted up in this way, attribute--value problems (AV problems~ o More precisely we may now introduce the following definitions and properties. Definition 3. I. In a problem, we identify a set A of elements called attributes: A - - { A I, A 2, -.., Ai,...,An~

(3.1)

C~

Definition 3.2. Each attribute A i is associated with a value set Vi, i.e. : i i •, vii~ V i = { v~, v2,...,vj,.. n where each v i is called value for the attribute A .

(3.2)

153

Definition 3.3. An attribute- value couple (AVC) for a parameter A i is a couple G i =(Ai, v~) where Ai is an attribute, ans v~ is a value for the attribute A .

iC]

Theorem 3. I.

An attribute A i individuates a set C i of AVC's, called AVC set

for Ai. Proof. 3y Definitions 3.1~ 3.2 and 3.3, we construct Ci in the following wa~v : i 0

We arrive now to the main definition in which the notion of structure of state is introduced.

a

Definition 3.4: A structured state (S S) N (i.e., a state with structure), is an n-tuple of AVC's, s.t.: = ~ c i, c 2 ,...,ci,...,c n ~

(3.4)

where each c i is an AVC f o r each parameter A i of A~ i.e. :

(VA i) ((Aifia)^(o i = (Ai, v~)))

(3.5) D

We may now introduce the notion of state space (structured). Theorem 3.2. The SS set (i.e., the set of structured states) S is the cartesian product of the AVC sets for all the attributes A of A, i.e. : 1

= C 1 x C2 x . . .

x C.z x . . .

x On

(3.6)

PRof. 9~2om t h e d e f i n i t i o n o f c a r t e s i a n p r o d u c t ) an e l e m e n t ~ o f S i s , of (3.6), an n-tuple~ made up with elements of each AVC set C..

because

l

From (3.3) we obtain therefore the SS ~, s.t. the (3.4) holds°

n

Please note that, within this description, we can clearly understand what does it means that two states s' and s" are different: it means that they have some AVC's (@.g.,c~ and c:'] which are different, i.e., some attribute A i l" takes two differente values (e°g., v i and v. ). J Now, in order to set up the notion of problem, we have to introduce somehow the notion of transformations from an SS ~ to an other SS s* (which in the syntacc tic description were formalised with the operators oflC). If any transformation between anJ two states would be possible, we would be faced by a very trivial situation in which any problem (i.e., any choice of initial and final states) would be solved with just an one-step solution. Problems of this n a t u r % called universal problems, would be represented by a complete graph, because all the transformations are possible, i.e. ~ any two vertices are connected by an arc. In this case, the set of the operators, in the syntactic description~ would be represented by the set of all the functions on S, i.e.:

where r u is the universal operator set. In a real problem~ on the other hand, we are faced with some restrictions with respect to the possible transformations~ therefore the graph associated with the problem is an incomplete graph, and the set of operators i~ would be less powerful than F u.

154 The idea which is embedded in the semantic description is practically this one: to describe in the formal way, the constraints which are imposed on the pos sible transformations between two SS's. Moore precisely, it will be exposed how to describe, with expressions which deal with attributes and values, i.e. ,with the structure of a state~ the limit a tions which make up a real problem (i.e. 9 an incomplete graph) starting from an universal problem (i.e., a complete graph). It is interesting to observe that these expressions will present two different mathematical aspects, namely, the algebraic one, and the logical one. More precisely we Definition 3. 5 .

may now introduce the following definitions and properties.

A legal condition ~LC)L i

L. ¢ I

or

is a binary relation on S, i.e. :

SxS

(3.8)

where Pi is a predicate, i.e. : Pi(s~,s") : S x S - (T a ~

,{T,F~

(3.10)

F stand for ~rue and ~ .

[]

Please note that in Definition 3.5, while (3.8) implies the algebraic nature of a legal condition, (3.9) and (3.10) illustrate its logical aspect: in particular the predicate Pi can be expressed within some logical calculus (e.g., the first-order predicate calculus) operating on the attributes and the values as variables and constants. A legal condition is therefore a mathematical expression, drawn from the intuitive notion of a problem, which enables one to describe some of the limitations existing on the transformations between SS's. Therefore a problem can be made up with a certain number of these expressions. Definition 3.6. An LC set L is the set : L ={%,

LI, L2, ..., Lh, ..., L j

(3.11)

where I L and

Lh, 1~h_~r

o

= S x S

(3.12)

are all the L0's.

When there are not LC's, we absume the existence of the special LC L o which yields the universal problem. With the next definition we introduce an important notion which takes, in the semantic description, a place equivalent to that of r in the syntatic description. Definition 3.7. The ognstraint O of a problem is a bynary relation on S (i.e., O C S x S), s.t. : r C = ~ L.,1 (3.13) h=O

[]

[]

The constraint C is therefore the u~ifioation of all the formalised information which makes up a real problem from an universal problem. This concept provides a new way of defining the notion of problem which

is

155

equivalent to that one given in the syntactic description. Definition 3.8. An AV problem schema M is a couple M = (S, C) where S is an SS set, and C is a constraint. C~ Definitipn 3.~. An AV problem P is a triple P = (MI i, f)(or P =(M, ~, K)) where is an AV problem schema, ~ is an initial SS and f is a final SS (or E is a set of final SS's). We can now illustrate the main result, which shows, essentially, the equivalence between the semantic description and the syntactic description. Theorem ~. 3. The set of operators F and the constraint C are equivalent. Proof. First of all we recall from algebra that~ given a function f on a set Q , i.e. : f , Q .----~Q ( 3. 14) f : q'~--~q" (3.15) we can associate f, with a binary relation G(f) on Q, called the graph of f,s.t.G(f) = I(q',q")i(q',q"EQ)A(q" = f(q')))

(3.16)

It exists clearly an equivalence between S and S, i.e., between s and N, because each SS ~ (i.e., each n-tuple of AVC's) constitutes the formal structure of each state s. Now, if we consider the set of operators F , because of Definition 2. I., and (2.1) we may obtain the constraint C, in %hisway : C = V

¥.61"

G(~i )

(3.17)

&

Please note that the C which is obtained from I~ is unique. On the other hand, if we consider the oonstraint C, which is a binary relation, we may set up a set of operators F in this ws~v :

In this case the choice of the operators y. of F consists in the covering of a binary relation with functions (i.e., the operators yi ) , which may be done in many different wa~vs, all equivalent to eaohother. ~S We may therefore state these two final results. Thcore m 3.4. An AV problem schema M is equivalent with a problem schema Mo Proof. From Theorem 3. 3 we obtain directly that B is equivalent to S, and C is e~Ivalent %0 rTherefore, from Definitions 2.1 and 3.8, we obtain that ~ is equivalent to M.

[] Theorem 3.5. An AV problem P is equivalent with a problem P. Proof. From Theorem 3.4 we derive that M is equivalent to M. Moreover, from Theorem 3. 3., since S is equivalent to S, we have also that ~ is equivalent to i, and f is equivalent to ~ (or K is equivalent %o K). Therefore, from Definitions 2.3 and B.9., ~we conclude that P is equivalent to P.

156

IV.

T~

C0~[PUTATION OF THE EVALUATION FUNCTION

We will briefly outline a method, based on the semantic description of SSPS, which prevides the computation of the evaluation function ~(~). This technique is based on the idea of computing the estimate ~(~), by solving an auxiliar,Z problem in which ~ is the initial state, and the solution is easy to be found, and provides a lower bound of the cost h(s) of the optimum solution for the main problem. Since all the auxiliary problems share the same goal problem, we shall intro_ duce the notion of auxiliary goal problem. Definition 4.1. Let F -- (S, C, 3) (or F = (~,C,K)) be an AV goal problem, obtair.n ed from an AV problem P',F'=(~', C', f~)(or F'=(S', C', KI))is called an auxiliar~ AV ~oal problem for F iff:

(i)

,~ -- ,~'

(ii)

f -- f'

(iii)

.T.,'.C T,

(4.1) (or K = K')

(4.2)

(4.3)

where L t and L are the LC sets from which C t and C are obtained. We shall indicate :

F' ~ ~"

(4.4) r-',

In other words, we say that F'~.F ~if they have the same SS set,the same initial and final SS's, and if the LC set L' is a subset of the LC set L. Theorem 4-I. If F ' ~ F then :

O'~O

(4.5)

Proof. The proof is obtained directly from Definition 3.7 and from (4.3) because of the proper~y of set intersection. Now we have the main result which provides the computation of the estimate. Theorem ~.2. If F ' ~ F then for all s ~ S , the optimum cost of the solution to the AV problem P~=(F',~)=(~',O',s,f') (or PI=(M',C',s,K')), is a lower bound ~(~) with respect to the AV problem P=(~,s,~)(or P=(~,~,~)). Proof. Because of Theorem 4. I, C'~ C 9 therefore the graph G=(S,A) and the graph G'=(S',A') respectively associated to the AV problem schemata N and ~' are such that, the set of arcs A is a subset of the set of arcs A' (see Definition 3.7 and Theorem 3. 3), i.e. :

~') A Then, if ~

(4.6)

is an optimum solution for the problem P , (i.e.,a minimum path

in a fro~ ; to ~ (or

~)), ~

is also ~ solution

~n

for the proble~

~' (i.

e., also a path in G', but not necessary a minimum one). Since h(s) is, by definition/ the optimum cost of a solution for P, then the optimum cost of a solution for P', which is not ne@essarily a solution for P, is a lower bound for h(s), i.e., it is the estimate h(s). Theorem 4.2 constitutes the basis for a new algorithm~ in place of the Hart 9 Nillson, and Raphael algorithm, for the heuristically guided search. The basic steps of the algorithm are the following ones: I. Given an AV problem P, and its associated AV goal problem F, an auxiliary AV goal problem F' is automatically constructed by eliminating some elements of the LC set L.

157

2. The computation of the estimate h(s) for the problem P is performed by solving the problem P'=(F'~,~). An evaluation of the complexity of this method, compared with the well known search strategies, will constitute the goal of a new research effort. V.

C 0 N C L U S I 0 N S

In this paper we have proposed a new description of SSPS, called syntactic description~ which constitutes a framework useful to structure a rich content of informations about a given problem. The semantic description has been shown to be equivalent to the syntactic description of SSPS. A method has been briefly outlined, which is based on the semantic description~ and which makes possible to extract~ in an automatic w~v, i.e.~ by computation, the heuristic information useful to guide the search in the state space. The future research work shall be addressed to the exploitation of the semantic description as a powerful basis on which to formulate procedures apted to solve the main goal of automatic problem-solving. In particular the direction of heuristic guided search and learning will be investigated. Acknowledgments The authors express their gratitude to Dr.D.Mandrioli~ for many stimulating discussions, and to all the researchers of the Milan Polytechnic Artificial Intelligence Project. REFERENCES Feigenbaum, E., aud Feldman, J.(eds). Computers and Thought. Mc Graw-Hill Book Company, New York, 1963. Minsky, M. (ed). Semantic Information Processin ~. The M.I.T. Press, Cambridge, Massachusetts, 1968. Nilsson~ N° Problem-Solving Methods in Artificial Intelligence. Mc Graw-Hill Book Company 9 New York, 1971. Slagle, J. Ar$ificial Intelligence: the Heuristic Programming A~roach.~'c GrawHill Book CompanY, New York~ 1971. Amarel, S. On Representations of Problems on Reasoning about Actions. Michie,D. led.) Machine I~telligence ]. pp. 131-171. American Elsevier Publishimg CompanY, inc., New York, 1968. Hart, P., Nilsson~ N. 9 and Raphael, B. A Formal Basis for the Heuristic Determination of Minimum Cost Paths. IEEE Trans. Sys.Sci. Cybernetics~ vol. SSC-4, no.2, pp. 100-107, July 1968. Mandrioli, D., Sangiovarnai Vincentelli, A., and Somalvico~ M. Towards a Theory of Problem-Solving. Marzollo, A. (ed). Topics on Artificial Intelligence. SpringerVerlag New York, Inc., New York, 1973. Sangiovanni Vincentelli~ A., and Bomalvico, I~{.Theoretical Formalisation of the State-Space Method for Automatic Problem Solving. (In Italian). MP-AI Project. MEMO $~-AIM-6, October 1972.

158

Sangiovanni Vincentelli, A.~ and Somalvico, Mo Theoretical Aspects of StateSpace Approach to Problem-Solving. Proco VII International Cor~ress on C~bernetics, Namur, September 1973 (a)o Sangiovanni Vincentelli~ A.~ and Somalvico M. Problem-Solving Methods in Computer-Aided Medical Diagnosis. Proc. XX International Scientific Conference on Electronics, Rome, March 1973 (b).

PERTURBATION THEORY AND THE STAT~.M~NT OF INVERSE P R O B L ~ S G. I •NARCHUK Computer Center, Novosibirsk 630090, U.S.S.R. Some aspects of the theory of inverse problems for linear and quasilinear equations of mathematical physics are investigated. The theory of perturbations is developed with respect to selected functio nals of problems. I. Conjugate Ftmqtions and the Notion of Importance Let us consider the function

(~.1) where ~

90(X) which satisfies the equation

k ~(x) = ~ (x), is some linear operator and ~

sources in the medium. Here I

(~',)

is the distribution

of the problem (time and space coordinates, us assume that the operator ~

of

is the totality of all the variables and functions

energy, velocity). Let 5o are real and

that

g~. For the sake of definiteness we assume, for example, that the pro cess under investigation is connected with diffusion or transfer of a substance though the conclusions of the theory are beyond the scope of this discussion. Let us introduce Hilbert space functions with the scalar product

where integration is carried out over the whole domain ~ of the functions ~ and ~ . The usual purpose of solving some physical problems is to obtain some quantity which is the functional of ~ (x). Any value linearly related to ~(az) can be represented as such a scalar product. For example, if we are interested in the result of the measurement of some process in the medium with the charactecteristics vice Z(,%') the value is

(~.3)

Jz =S~(x) 2

of the de-

( x j d ~ = (~t, z ) .

Thus we consider the physical quantities which can be represented as a linear functional of ~ (X)

,7e £ ~ J = (~, , ) ,

160

where the quantity p

is a characteristic of the physical process

under consideration. Let us introduce along with the operator ~ its conjugate operator ~* , defined by the Lagrangian identity

(1.4)

(~7,

LAJ

CA,/~#)

=

for any functions ~ and ~ Together with (1.1), which will be called the basic equation, we introduce, first formally, a conjugate inhomogeneous equation

where p (~c) is some arbitrary function and ~@ ~ ~ ~. Substituting solutions (1.1) and (1.5), ~ and ~ into (1.4) in place of the functions ~ and ~ , we obtain

(n.e)

(~p,

=

D ~p )

or, using equations (1.1) and (1.5), (1.7)

~'cpJ, ~,; = c~,, °-~,

i.e.

Therefore, if we want to find the value of the functional ~ [~] we can get it in two ways : either by solving equation (1.1) and determining this value according to the formula

(1.8)

Jp [ ~]

=

(~, p ,

or by solving equation (1.5) and determining the same value according to the formula

(1.9)

Jr c ~7 = J~ [ p]-- ~p, ~j.

Consequently, the function ~p (~) , satisfying equation (1.5) can be put to correspond with each linear functional JpL'CP.7= (9~,.o ) where the function p¢~C) , which characterizes the measuring device, should be used as the free term of this equation. 2. Perturbation Theory for Linear Functionals If the properties of the medium with which the field interacts change i.e. if the operator of equation (1.1) becomes

161

then both the field

~C~J

and the value of

Jp E W ]

change :

Let us find the relation between the variation of the operator ~ L and that of the functional SJp . The perturbed system is described by the equation

(2.1)

L~'-- (z.+ 5D)~,~= ~'

The conjugate function of the unperturbed system corresponding to the functional ~p is described by the equation

L * ~p~

(2.2)

= p.

Multiplying scalarly equation (2.1) by ~ , equation (2.2) by ~i and subtracting one from the other, using the determination of the conjugate operator of (1.4) we obtain, on the left,

(2.3)

c~,~ z,'~,'j - c~,~,z

(v3)

Input

vJ I

HYDROLOGICAL

Noise

SYSTEM

( vl , v2 , ~(in) ) i

ERROR

I Measurement C~m)

CRITERION

u

~0(k+l) = ~(vl(k) , ~0(k) ) ~l(k ) = ~i ( £0(k ) )

1

I Simulated L measurement

MODEL

ADJUSTMENT

£2 = ~2 ( ~ ' ~i (k)) A~

S(out) = ~(~,S(in),~l(k),~2,~)

(O(OuT))~°utP ut

Fi~. 8.

SimulaYion and adjustment plan.

STI~TE~Y

Y

o

60

150

,.L.i,~l,,.,,k

120

i

J

180

rJ

j

210

/

210

~'

I r.,,.. I, \


Z

/

/

I~u

V

Fig.

5

Time evolution for the staggered

scheme.

Moreover for this scheme a set of equations) equation (7)) have been used:

u

m,1

he

-G-E-* . ~

equivalent

to

v'.t '/:

v {u'+ v"3 k~_

(19)

: o

In this set of equations (19) u) v are the x and y mass flux. The advective terms are not taken into account because of their small influence~ especially with the present network.

210

Moreover the elimination of these advective terms allows the construction of a second oPder scheme (cfP. (6), chapteP III, paPt i). However we must take care in the computation of the last term which takes into account the loss of enePgy due to the friction (7). This term must be evaluated durin~ the second step as descPibed at the beginning of this paPagPaph, but at the time t+At only the water levels are known. For that Peason we must wPite the finite difference equation in the following form: y~+2

u,j .......

z

j

The boundaPy conditions which define the vanishing of the normal component of the velocity to the close boundary, ape imposed only at these time levels:

The water levels on the open boundaries are imposed at the time levels:

.4)

About some simulations of the pPopagatio n of the tidal way e inside the Venetian Lag0o ~.

FoP this method too the stability analysis in the fPequency domain was performed. The scheme is stable even if the period of the input wave is of 1 hour. The simulation of the semi-diurnal tide showed at the point 6 (fig. 4c) a gain factor (Patio between the tide inside and outside the lagoon) of 1.2 with a good agreement with the expePimental data. A few remaPks about the convergence of the numerical scheme. Four 9-day simulations have been carried out with a 1250m-size mesh network and with time step At 10,20,30, $0 seconds. In oPdeP to illustrate the rate of convergence of the numerical solution n(At,x,y) to the true solution ,(x,y,t) we have plotted in fig. 6 the function F(t,At) which is a sort of measure of the space averaged deviation of ,(At~x,y,t) fPom ~(lO,x,y,t)

D is the whole sumface of the lagoon ,(lO,x,y,t) is the computed water level with a time step of i0 sees n(dt;x~y,t) is the computed water level The graph in fig.7 clearly shows the following pPoperties: a) the magnitude of the errors is dependent on the Pate of vaPiation of the water level but it seems to be limited b) the maximum mean erPor is of 3 m m

211

.5

.0S

24,

:/2

Fig.

6

An example of the convergence

96 of the numerical

,lo~

solution.

9

Fig.

7

A 9-day simulation with the 1250m-size mesh network.

An example of simulation is shown in fig.7 (Ax = Ay = 1250m , At = 60 see). The roughness coefficient c was set to 60 over the whole lagoon. The simulation shows a good qualitative agreement with the experimental data. For instance the oscillation in the northern part of the lagoon (lines 8, 9) is strongly damped as in the field. From the inspeotion of the right part of fig.8 we can evaluate the phase of the tidal wave at different points of the lagoon: the agreement with the experimental data is satisfactory as the relative error never exceeds 20 percent. Another simulation was carried out with a network of 1600 nodes (mesh size = 625 m). The main results of it were an encouraging decrease of the relative error in the phase evaluation and a better d e s c r i p t i o n of the tidal wave in the northern part of the lagoon. This fully justifies our efforts to build a new network (300 m~ for instance) and to analyze a sufficien~amount of mareographic data in order to obtain a reliable calibration of the model.

212

References I.

J.J. Lee, Wave induced oscillations California Institute of Technology, Dec. 1969.

in harbors of arbitrary shape, Pasadena, Report No. KH-R-20,

2.

J.J. Leendertse, Aspects of a computational model for long-period water-wave propagation, Rand Corporation, Santa Monica, Ca., RM5294-PR, May 1967.

3.

J.J. Stoker, The formation of breakers and bores, Comm. Pure AppI. Math., i,I, 1948.

4.

G.F. Carrier, H.P. Greenspan, Water waves of finite amplitude on a sloping beach, J. Fluid Mech., 4,97, 1958.

5.

P. Lax, Weak solutions of non linear hyperbolic equations and their numerical computation, Comm. Pure Appl. Math., 7,159,193, 1954.

6.

S.K. Godunof, V.S. Ryabenki, The theory of difference schemes North-Holland Publ. Co., Amsterdam.

7.

R.D. Reid, Numerical model fop storm surges in Galveston bay, J. of the Waterways and Harbors Division, February 1968.

SEA

LEVEL

PREDICTION

by A. Artegiani, CNR ^ IBM

MODELS

A. Giommoni,

FOR

VENICE

A. Goldmann,

P. Sguazzero,

Laboratory for the Study of the Dynamics Scientific Center, Venice

A. Tomasin

of Large

Masses,

:~ Venice

In the northern side of the Adriatic sea (Fig. I) the water level, under certain weather conditions, can rise beyond the ordinary value of the high tide and cause severe floods in the coastal areas. The effects of these phenomena are particulary famous for Venice, but they threaten a much larger area, which certainly includes the whole lagoon where Venice is located. Venezia

~f~vinj

Marina di~ A a v e f ~

l~ul~

Ancon \.

"X e eark

~ ,

Split

.Subrovnik "X,

~i~anoA~ used in the present work.

The possibility of forecasting the floods was limited, until recently, to the experience gained by the local offices in the past, and enabled to a prediction m a d e two or three hours in advance.

ea

214

But an objective forecast with a longer time lag is necessary, and the need will increase when sluices for the lagoon, which have already been planned, will be built. From a practical point of view, we can observe that winter time is the period in which the danger for high water in Venice becomes most threatening. It is in this period that cyclones coming from the British Isles or from the Western Mediterranean (Fig. 2), move towards Italy.

- F i g . 2 P r e f e r e n t i a l t r a c k s of A t l a n t i c c y c l o n e s e n t e r i n g the M e d i t e r r a n e a n

These atmospheric disturbances generate winds blowing from the South or the South-East which pile up water in the northern part of the Adriatic Sea. An experienced weatherman, just looking at the meteorological maps, can guess whether dangerous storm surges will occur or not, but if we want to give a numerical forecast of the sea level with the explicit indication of the time of the event, a more careful analysis is necessary. The observed level at any point in the Adriatic can be regarded, in a first approximation, as the sum of two components which we assume to be indipendent: a) the astronomical tide due to the attractions of the sun and the m o o n on the waters of the sea; b) the meteorological tide due to the atmospheric forces, n a m e l y the tangential wind stress and the n o r m a l barometric pressure. A s the astronomical tide can be predicted with high accuracy with standard techniques w e shall speak f r o m n o w on only about meteorological tide.

215

Given the geometric characteristics of the Adriatic, the surges in the sea can be described by the h y d r o d y n a m i c equations in one-dimensional form:the continuity equation which follows f r o m the principle of conservation of m a s s

?--g *

ax

the m o t i o n e q u a t i o n w h i c h s i m p l y e x p r e s s e s the b a l a n c e b e t w e e n the t e r m s of inertia and surface slope and the atmospheric pressure gradient, the wind stress, the bottom stress

where X

is t h e basin

axis

Q [×,~) the discharge ~ [X,t~ [X~ ~-[x~ k[ x) ~ ~'$ ~

the displacement of water from equilibrium surface the cross section area the basin width the basin depth from equilibrium surface the atmospheric pressure the surface wind stress the density of the water the acceleration of gravity T h e bottom stress T h is usually simulated by a quadratic or linear function of the m e a n velocity of the Yluid: given the use w e will m a k e of this kind of analysis the choice of a frictional t e r m proportional to the discharge is satisfactory. It is easy to see that in this case the solution of the h y d r o d y n a m i c equations (under proper initial and boundary conditions) takes the f o r m of a series of fundamental oscillations (modes): each of t h e m can be obtained as a convolution integral of a suitable forcing function (depending on the me[eorological variables) and a response function of a sinusoidal-damped type. The behaviour of the Adriatic Sea s e e m s sufficiently close to this simplified scheme; it's k n o w n f r o m m a ny years that the basin oscillates on characteristic frequencies, when it has been excited by atmospheric disturbances. The first eigenfrequencies correspond to the dominant oscillations or seiches, whose shapes, derived from the previous theoretical considerations, are shown in figures 3 and 4. The first shows the uninodal seiche, with the node at Otranto, the mouth of the basin, and with a period of about 22 hrs. The second figure shows the binodal seiche which as a period of about I0 hrs.

216

/..-r

Elevation I 0~ (arbitrary- uni~sll

/

f T

~

I

Ill

;

;.

Distance (Km)

- Fig. 3 Theoretical profile of the first seiche along the Adriatic. Experimental estimations are also given I00 50

-50

~

~

~ ~

"~ "~

o

~o

~ ;>

,.Q Q

Distance (Km)

- Fig. 4 Theoretical profile of the second seiehe along the Adriatic Experimental estimations are also given

Spectral techniques were used in order to obtain an experimental evaluation of the periods and put in evidence the importance of the seiehes. The spectrum of many months of the recorded sea level at Venice, 1966, from which the astronomical tide has been subtracted is shown in Fig. 5. Apart from a very low frequency band which is due to long period meteorological fluctuations, the spectrum contains two distinct energy bands (]3, C) at the seiches periods of 21.7 and I0. 8 hrs; the periods computed with the simple ana lyrical model previously presented are of 21.69 and 9. 92 hrs, and the qualitative agreement is good indeed.

217

10 2 Energy Density (cm2/cph)

I0 0

10 - 2

DB

d. 03 -

t ~

S I

o" 06

C ,I

o. 09

,

0:12

o.

p

Is

Frequency

(cph)

A more direct analysis can be carried out by filtering the time series sea level of several harbours in the Adriatic Sea by suitable band-pass ters centered on the eigenfrequencies and by comparing the outputs.

of the fil-

Fig.

5 Spectrum

of t h e m e t e o r o l o g i c a l

Elevation

tide in Venice,

1966

(cm)

Venezia

Rovinj -20 Ancona

201 - 20 :

"~

Ortona -20 20 Split

-20

]Vianfredonia

Dubrovnik

Otranto

201

-2(? 201 -20:

r~

20 -20 0

24 48

72

96 120 144 168 192 216 Time (hours)

- Fig. 6 Simultaneous Progressive hours

plots of the seiehe in different from 16 February 1967, 7 a.m.

harbours

of the Adriatic•

218

Ing Fig. 6 are shown the plots of the first seiche along the Adriatic. According to the theoretical deductions, the amplitude of the seiche is decrea sing from the North to the South; furthermore it appears evident that the seiche is damped in time as an effect of the friction. The most important assumption made in the previous exposition was the linea rity of the surge; this means that the equations governing the motion are linear with time-indipendent coefficients. It follows that the solution can be obtained by applying a suitable integral operator to the forcing functions namely the wind stress and the pressure gradient:

The p r e d i c t i o n p r o b l e m for the s e a l e v e l ai the point x 0 (say Venice) is to find the function R (x~, x, T) (the r e s p o n s e function) w i t h r e l a t e s to ° ~[.CX.,'i') i t s principal causes :~ (x,T) W e approach the problem of determining the response function of the Adriatic, with statistical methods. Essentially the following steps are involved: a) the integral operator is discretized

b)

the prediction weights R., are determined & la least squares from discrete .IK time series of input (drlwng functions) and output (recorded sea level at x 0, (say Venice again) corrected for the astronomical tide). Let ud recall that the problem we are interested in is to furnish a short time prediction of the sea level at Venice, so that the formula in which the prediction weights, are to be determined becomes

where d is the forecast interval and the F(J)'s are the actual predictors. The definition of the forecast interval, set to 6 hrs, was suggested by experience and statistical analyses. The cross correlation curve between the meteorological tide in Venice and the estimated wind stress in the middle Adriatic (Fig. 7) indicates a mean delay of 6 hrs between the two time series. As far as the choice of the input variables is concerned, we decided to use the following groups of predictors: I) values of pressure gradient in the Adriatic interpolated from coastal pres sure values of 7 meteorological stations, namely Venezia, Marina di Ravenna, Bari, Pula, Split, Dubrovnik (see Fig. I); for this purpose the whole basin was subdivided in 5 parts in which a mean gradient was assumed;

219

Cross correlation coefficients

0.4 0.3 0.2 0.1 0 ,.

. . . . . . . . . .

-0.1 -0.2 12

24

36

48

60

72

, 84

, 108

96

• * 120 132 Time lag (hours)

- Fig. 7 Cross meteorological

correlation functions for the wind tide in Venice (1966)

stress

in the Adriatic

and the

2)

values of IW~VP for the same zones; they are representative of the surface wind stress by accepting the geostrophic hypothesis and a qua dratic relationship between stress and wind. It can be shown that for the meteorological tide the autocovariance remains large for a considerable amount of time. In other words the meteor_o logical tide is highly predictable from its own past. Se we decided to introduce as an additional input in the set of predictors 3) the past of the sea level itself 80 Elevation

(cm) 60

40 20 0 -20 -40 12

24

36

48

60

72

84 Time (hours)

- Fig. 8 Observed (~) and predicted (---) meteorological Progressive h o u r s f r o m 1 A p r i l 1971, 7 p . m .

tide in V e n i c e .

220 Using 4 years of mareographic and meteorological data, we obtained the follo_ wing results: we were able to account for 87.5 percent of the variance of the meteorological tide, leaving a standard error of 5.7 cm (in terms of the observed sea level, including the astronomical tide the predicted variance was 96.5 percent). The model was applied during all cases of high water, starting from autumn 1970. In figures 8 and 9 two cases of high water are shown, in which we have plotted only the meteorological tide, because the computation of the astronomical tide is standard: the solid line represents the surge as it can he observed, the dashed line the surge forecasted 6 hrs in advance using the statistical model. 8O Elevation

(cm)

t

40

i

20

-20 t V

. . . .

-40

, 12

24

36

48

60

72 Time

-Fig. 9 Observed (--) and predicted (---) meteorological Progressive hours from 12 February 1972, 4a. m.

(hours)

tide in Venice.

Following more or less" the same line, another approach was developed, which directly performs the integration of the basic equations (water conservation and motion) that were shown above. Since the description of the meteorological "forces" is available, as it has been shown above, and the geometrical characteristics of the sea are obviously known, a simple numerical technique was applied, using a computer. An explicit finite-difference scheme took the place of the differential equations, and it was clearly a one-dimensional model, since, as stated above, the Adriatic can be considered as a channel, with varying depth and cross-section. (For special purposes, also a two-dimensional model was used, but for the practical task of surge prediction the former one was sufficient). When using this method, one can integrate the equations and deduce the levels and currents only as long as the forcing functions are known. To predict future conditions, one has either to use the weather forecastings or just to keep the latest weather reports as valid in the following hours.

221

This surprising "wind freezing" turns out to be useful when six hours in advance for the reason that this lag corresponds response time (in the obvious sense), so that the changes in cal field that will occur during this time will not appreciably is. This method also was tested by a "field use", and again the tisfactory, as one cas see from Fig. I0, whose accuracy is

forecasting about to the typical the meteorolo~ affect the resul results were typical.

sa-

120

Elevation

(cm)

J

110 100 90 80 70 60 50 40 30 20

...........

0

3

6

9

12

15

18

21

24 Time (hours)

numeri - Fig. 10 An example of prediction of floods in Venice by hydrodynamical cal models. Using the meteorological reports of the time marked on the horizon tal axis, a forecast was issued (square points) announcing the flood, as soon as the reports were made available (more than an hour later). Both the solid line (observed values) and the prediction refer to the 'total' tide, without astronomical subtraction. The time scale begins at 9. 0, November 9, 1971 The conclusion that one can draw from our experiment is that the forecast of the floods in the Adriatic can be made, with a good accuracy, up to six hours in advance. A longer time lag for prediction will require a deeper involvment in the meteorological forecasting. We are indebted to the Ufficio Idrografico del Magistrato alle Acque di Venezia, the Servizio Meteorologico dell'Aeronautica Militare Italiana, the UNESCO and many of our colleagues for their help in this research.

OPTIMAL

ESTUARY

AERATION:

AN

A P P L I C A T I O N O F DISTRIBUTED P A R A M E T E R CONTROL THEORY Wayne Hullett* The S i n g e r C o m p a n y Librascope Systems Division ABSTRACT T h e concentration of dissolved oxygen in a river has c o m e to be accepted as a criterion of water quality. In-water tests have s h o w n that artificial aeration by m e a n s of in-stream mechanical or diffuser type aerators can be an economically attractive supplement or alternative to advanced wastewater treatment as a m e a n s of improving water quality. This paper applies distributed p a r a m e t e r control theory to obtain the aeration rate that m a x i m i z e s the dissolved oxygen distribution with the least control effort. Both the s y s t e m state and the control input are distributed in space and time. A m e a n square criterion functional is used which allows the optimal feedback control to be determined as a linear function of the state. T h e feedback gain is found as the solution to the infinite dimensional equivalent to the matrix Riccati equation. A n analytic solution for the feedback gain is found for the non-tidal portion of the river, which is modelled by a first order hyperbolic equation. T h e estuarine portion is described by a diffusion equation, and a n u m e r i c solution obtained by approximating the diffusion equation with a finite dimensional system. A n e x a m p l e is given using historical data f r o m the D e l a w a r e estuary, and the dollar cost of the optimal control is c o m p a r e d with other ad hoc control strategies.

1.

INTRODUCTION

T h e concentration of dissolved oxygen (DO) in an estuary has c o m e to be accepted as a criterion of water quality.

It has been suggested (Whipple,

1970) that the use

of artificial in-stream aeration can be an economically attractive supplement or alternate to advanced wastewater treatment. T h e p r o b l e m to be addressed here is the determination of the aeration rate that achieves the best increase in the D O

level with the least control effort. This is

approached as a distributed p a r a m e t e r control problem,

w h e r e both the s y s t e m

state and the control are distributed in space and time.

A solution in the f o r m of a

feedback control is sought. Z.

WATER

QUALITY

MODEL

T h e dissolved oxygen level in an estuary is decreased by the oxidizing action of an effluent, described by its biochemical oxygen d e m a n d (BOD).

The DO

by atmospheric aeration and by artificial in-stream aeration.

Other factors affect-

ing D O

but not included here are benthal demand,

and the nitrogen cycle.

algal respiration,

is increased

photosynthesis

In order to illustrate the action of the o p t i m a l control,

the following partial differential equation (Harlemar~

1971) is used to describe

* T h e research described in this paper w a s conducted under the supervision of Professor A. V. Balakrishnan while the author w a s a R e s e a r c h Assistant at the University of California, L o s Angeles.

223

t h e d i s t r i b u t i o n of d i s s o l v e d o x y g e n : bC

1

b

bC

5C

~x = ~ ~x (AE ~ x ) ' V ~ x

(I)

+K3(cs-c) + u - ~i L

with initial condition C ( x , 0)

=

Co(X)

and b o u n d a r y c o n d i t i o n s

c ( 0 , t)

-- f(t)

C(xf, t) = g(t) where A(x,t)

= estuary cross

C(x,t)

= dissolved oxygen concentration(rag/C)

Cs(X, t) E(x, t)

s e c t i o n a r e a (sq. f t . )

= DO s a t u r a t i o n v a l u e ( r a g / g ) = t i d a l d i s p e r s i o n (sq. ft. / d a y )

IK1 = BOD o x i d a t i o n r a t e ( 1 / d a y ) K 3 = atmospheric

aeration rate (1/day)

L(x, t)

= biochemical

oxygen demand (rag/g)

U(x,t)

= aeration rate (rag/g/day).

T h e BOD c o n c e n t r a t i o n a s a f u n c t i o n of the e f f l u e n t i n p u t o b e y s a s i m i l a r bL i 3 ~L ~ " = X 7x ( A E T x ) -

?L v-~'x - K I L +

equation:

J

where J(x, t) : effluent discharge rate (rag/%/day M a n y criterion functionals are relevant: the least squares has the advantage that the optimal control is a feedback control.

T h e functional to be m i n i m i z e d is

therefore:

S(U)

=]

/

Q(x) [ C - C s] d x d t + X

J0 J0

R(x) U 2 dx dt J0 Ja

w h e r e t h e c o n t r o l i n t e r v a l (a, b) is c o n t a i n e d i n t h e e s t u a r y l e n g t h (0, xf), and T is t h e t i m e d u r a t i o n of t h e c o n t r o l . M i n i m i z i n g t h i s f u n c t i o n a l o v e r U(x, t) w i l l t e n d to d r i v e t h e DO l e v e l t o w a r d the saturation value, while penaiizing large control inputs. n e g a t i v e ) and R(x) ( p o s i t i v e d e f i n i t e ) a r e a p p r o p r i a t e

T h e f u n c t i o n s Q(x) ( n o n -

weighting functions.

224

3.

OPTIMIZATION

PROBLEM

Let

C(x, t) ~ L z (0, xf; 0, T)

= H 1

U(x,t) ¢ L 2 (a,b; 0, T)

= H 2

and define the m a p p i n g B:

H2--~H 1

by

I U(x,t);

xe (a,b)

BU 0;

otherwise

T h e m o d e l e q u a t i o n (1) c a n b e w r i t t e n i n t h e f o r m

~C --~-[- = F C

+BU+V

where F is the operator

1

~

=

X ~

with domain

D(F)

FC

(AE =

d e f i n e d by

~C ) ~C ~x -v~'x " K3C

[CCHllFC

CH1;

C(0, t) = f(t);

C(xf, t)

= g(t) ]

and V

= K3C s -K1L.

The criterion

functional is the sum of inner products in H 1 and H 2

The feedback control that minimizes u

1

=

1

- ~[ a

B*

S(U) ( L i o n s , 1968) is

(MC+G)

w h e r e M is a self adjoint linear operator and G(x, t) is a function that satisfies .~-

= -F,',~ M - F M -

~G

=

M B R -1 B~~ M + Q

and -F*

G +

:~ M B R -1 B;:,G - M V + Q C s .

T h e e q u a t i o n for M is the infinite d i m e n s i o n a l analog of the m a t r i x where equality M 1 = M 2 is to be interpreted M1C

= M2C

for allC

¢ D(F)

in the sense

(2) Riccati equation

225

T h e final conditions are MC(x,t)l

= G(x,T)

= 0

! t=T and boundary conditions are

MC(x,t) I

= MC(x,t) I

i x=0

= G(O,t) = G(xf, t) = O.

Ix=xf

Since final conditions are given for the equations for M and G, these m u s t be solved b a c k w a r d s in time,

so that knowledge of all future values of the s y s t e m p a r a m e t e r s

A, E, v and C s and the s y s t e m input V = K 3 C S - K I L This is not as bad as it seems,

however,

for t ¢ (0, T) is required.

and for the systems under consideration,

is actually a benefit, because it permits the use of all awailable information in computing the optimal control.

In general,

one is able to predict several days in

advance f r o m meteorlogic conditions and production schedules the temperature and BOD

discharges upon which C s and L depend.

it is reasonable that a control based on m o r e

Since this information is available, information would be better than one

c o m p u t e d without this extra knowledge. O n e can envision a system,

then, that would c o m p u t e an optimal control ( M and G)

for a period of, say, I0 days, for which the temperatures and discharges could be reasonably predicted.

When more

data b e c o m e s

available, after 5 days for

instance, a n e w M and G could be c o m p u t e d for the next i0 days, or for the remaining 5 days if it is desired to stop control {aeration) at that time.

O n e then has, at

each time, a control s y s t e m that is optimal based on all information available at that time. 4.

ANALYTIC

SOLUTION

T h e differential equations for M and G are extremely difficult to solve, and usually only a n u m e r i c solution is possible.

In the special case of a s t r e a m (no tidal action)

with a constant cross-sectional area and velocity, and a s s u m i n g the molecular dispersion to be negligible c o m p a r e d with the bulk transport due to the s t r e a m velocity, the differential operator F b e c o m e s FC

= -v

~C ~x

- K3C

In this simple case, an analytic solution for the feedback gain and forcing function can be found (Koppel and Shih, 1968). MC

= P(x,t)

C(x,t),

and for definiteness letting Q(x) = I; a =0; T

= xf/v

R(x) = I b =xf

Assuming

a solution of the f o r m

226

The Riccati equation becomes 8P ~P ~ pZ_ -~- = -v ~x + Z K 3 P + 1 with P(x,T)

= P(xf, t) = O,

which has the solution

I k[atanh(-(zt

P(x, t)

+~) - 143 ] ;

xvt

where

4/-K-K +7[ 2

=

I

= a T +tanh-1 ( ~ ) "~ = ve~ xf+tanh-i ( K3)~ " Substituting this solution into (Z) and solving (noting that the subsidiary conditions are G(x, T) = G(xft) = 0) yields T

i{A sinh(e~IT- o-] ) V

seth (-at+~5)

- cosh(-a0-+~)Cs[X-v(t-0-),

G(x, t)

-

]

~]} do-;

x< vt

=

vl sech(-~x+-?)

f xf{ o~ 4~sinh [v(Xf-O- ) ] V[°-'t-x'°'v ]

x

- u ] } du; _cosh(_ ~ ¢+.D Cs [~, t _ x--J"

x > vt

These solutions for P(x, t) and G(x, t) completely define the optimal control as a feedback function of the current D O distribution: u ( ~ , t)

1

= - ~ ( P ( x , t) C ( x , t) + G ( x , t))

227

5. In order to obtain results was approximated

NUMERIC: APPROACH

for the general

by a finite system

estuarine

case,

the diffusion equation

of d i f f e r e n c e e q u a t i o n s ( S a g e , 1968) a n d t h e d i s -

crete optimal feedback control found in the usual manner. This approach was applied to an example using historic Estuary

f o r t h e y e a r 1964.

The estuary

was partitioned

feet each and the control was constrained a b o u t 38 m i l e s

of estuary.

in sections

data from the Delaware i n t o Z3 s e c t i o n s

6 t h r o u g h 15,

T h e w e i g h t i n g f u n c t i o n Q(x) i n t h e c r i t e r i o n

w a s a d j u s t e d t o p e n a l i z e DO d e v i a t i o n s f r o m t h e s a t u r a t i o n

of 20, 000

representing functional

value over the same

interval. T h e o p t i m a l c o n t r o l f o r a 15 d a y p e r i o d i s s h o w n f o r s e v e r a l and the resulting

DO l e v e l s i n F i g u r e

control in Figure

3. 6.

The energy required

Energy

in horsepower

I - LTt - b = J o j a U(x,t)A(x,t) *7 (x, t)

d x dt

where

values (Whipple,

and pressure

in anaerobic

d e p t h of t h e a e r a t o r . was used.

H--p--Hr]"

1970) o f ~ 7 ( x , t ) f o r s t a n d a r d

c o n d i t i o n s of t e m p e r a t u r e

water vary from 0.68 to 1.36

LbO 2 depending on the Hp-Hr For the depths under consideration, a n a v e r a g e v a l u e of 1 . 0

This value can be converted

P zg. gz (Cs)T- C ~ =

1,

is given by

N (x, t) = oxygen transfer efficiency Measured

sections in Figure

w i t h t h e DO l e v e l s w i t h o u t

COST OF AERATION

for aeration

hours,

Z c a n be c o m p a r e d

to test conditions by the relation

T -Z0 c

(Cs)zo

r~T where P T

C

= pressure

( i n c h e s of m e r c u r y )

= temperature = empirical

(Deg.

C)

c o n s t a n t (1. 0Z5)

(Cs)20

= DO s a t u r a t i o n

value at standard

(Cs) T

= DO s a t u r a t i o n

value at test conditions

C(x, t)

= DO c o n c e n t r a t i o n

at aerator

T h e c o s t of t h e o p t i m a l c o n t r o l e x a m p l e other control schemes

i n T a b l e 1.

c o n d i t i o n s (9. 0Z m g / ~ )

f o r 15 d a y s a e r a t i o n

is compared

with two

228

1.0

o.9 \ ~ . , P. 0.6 ¢

0.5

w

0.2

x

0.1

"SECTION

SEC',ON,

_,

1

"" "- ~"~---~

1

.

0 0

1

2

3

4

5

6

7

8

9

10 1 1

12 13

14 15

TIME (DAYS)

FIGURE I .

OPTIMAL C O N T R O L (k=10)

9,0 SECTION 1 8.0 7.0 ~ Z '"

6.0

×

5.0

~ . . ~

....-,~ "" .L

O

"*'L ".,,,.,.. ~~ ..--- ~ , < : ~ ~ _ =..,¢ ~ " - ~ !

~

_.,,~ ~F ..-,.~'~,..~""~'~ 5 , ' ~....r

a

:::==1 SECTION SECTION ~ SECTION SECTION SECTION

i ~...--k"'1~

(3 4.0 w ->~ 3.0

~

.--- _ ----..

~

~

--~ ~ ' : "~ . . . .

~"

I l

--% ~'''~

5 7 9 19 17

SECTION 11 SECTION 15 SECTION 1•3

2,0

0

1

2

3

FI GURE 2.

4

5

6 7 8 9 TIME (DAYS)

t0

11

12

13 14 15

RESPONSE TO OPTIMAL C O N T R O L (k=10)

9.0 SECTION I

"3 8.o L

--~ 7,0

= I

~

~

SECTION 5

r......~______.~_.~.~__~.~b.~._~__..~...~.....,__.__ -SECTION 7

Z 6.0 UJ ¢3 ~ ~ > ~ X 5.0 _ ~ .._..._ ~ O

""

~-"~L-..,.--,-.~==',¢l

W~- ..m=~

iI~ ~

~ ~.,.~~ ..._.

~..__

4.0 ,

SECTION 19 SECTION 9 SECT1ON 17

w

>

3.0

~'~"",-..-~,-'--"

z.o O

~ ,-....~ ~

~,..~__.~,..~.~,.. ,.,__,.....

1.0

~

,,....-.-~"'~"""~

~,,......,~,~.....--,---.. /

0

1

2

3

4

FIGURE 3.

5

6 7 8 9 TIME (DAYS)

10 11

UNCONTROLLED DO

12 13 14 15

SECTION 15 SECTION 11 SECTION 13

229

COMPARISON

TABLE 1 OF CONTROL

STRATEGIES

Criterion Functional

Average DO Level (nag/C)

E n e r g y Cost (Dollar s)

Savings (Dollars)

Optimal Control

O. 115

5. 5

89,327

--

Control I

O. 168

5. 5

145,351

56,024

Control Z

O. 156

5. 5

105,046

15,719

Control 1 is an aeration rate that is constant over the s a m e distance and duration as the optimal control.

T h e constant rate is 0.5 m g / C /day.

Control Z is an aeration

rate that is proportional to the initial D O deficit in each section and is constant over the fifteen day duration.

The proportionality factor is 0.116.

T h e criteria used for c o m p a r i n g the effectiveness of the control strategies are the value of the criterion functional and the average D O trolled sections and the time interval.

concentration over the con-

T h e cost for electrical energy w a s esti-

m a t e d at $0. 006 per kilowatt hour, which is an approximate rate for large industrial users, it is seen for this e x a m p l e that the savings in energy costs obtained by using optimal control over the other controls is about $i000 per day. may

the

Of course there

exist other control strategies for which the savings would be less, but the two

chosen for c o m p a r i s o n are reasonable.

In this e x a m p l e the optimal control

approach is certainly a feasible alternative. 7.

CONCLUSIONS

A potential savings in energy costs of estuary aeration using an optimal control approach has been demonstrated.

T h e control of the aeration input rate can be

i m p l e m e n t e d as a feedback control, which has well k n o w n advantages in t e r m s of the sensitivity to incomplete knowledge of or variations in system p a r a m e t e r s

or

input s. Future plans include c o m p a r i s o n s with other sub-optimal control schemes,

other

n u m e r i c approaches and application of this approach to the control of other estuary variables such as the location and timing of B O D

discharges.

230

8.

REFERENCES

H a r l e m a n , D. R . F . "One D i m e n s i o n a l M o d e l s " , E s t u a r i n e Modeling, an A s s e s s m e n t , U . S . G o v e r n m e n t P r i n t i n g Office, Washington, D . C . , 1971, W a t e r P o l l u t i o n C o n t r o l R e s e a r c h S e r i e s , 16070 DZU 0~/71. Koppell L . B . and Shih, Y . P . " O p t i m a l C o n t r o l of a C l a s s of D i s t r i b u t e d P a r a m e t e r S y s t e m s with D i s t r i b u t e d C o n t r o l s , " I & E C F u n d a m e n t a l s , Vol. 7, No. 3, pp. 414-4ZZ, 1968. L i o n s , J . L . C o n t r o l e O p t i m a l De S y s t e m s G o u r v e n e s P a r D e s E q u a t i o n s A u x D e r i v e e s P a r t i e l l e s , Dunod and G a u t h i e r - V i l l a r s , P a r i s , 1968. Sage, A . P . O p t i m u m S y s t e m s C o n t r o l , P r e n t i c e Hall, Englewood Cliffs, New J e r s e y , 1968. Whipple, W . , et al " O x y g e n R e g e n e r a t i o n of P o l l u t e d R i v e r s , " E n v i r o n m e n t a l P r o t e c t i o n Agency, Washington, D . C . , 16080 DUP 12/70, D e c e m b e r 1970.

INTERACTIVE

SIMULATION FLOOD

PROGRAM

ROUTING

FOR

WATER

SYSTEMS

F . G r e c o a n d L. P a n a t t o n i IBM Pisa Scientific Center

- Italy

Introduction A computing system channels,

is here presented

with special regard

to flood routing in rivers.

built up during the study for the construction river basin.

It a l l o w s u s t o r e p r o d u c e

artificial watercourses.

f o r t h e s t u d y of w a t e r f l o w i n o p e n

The system

This system

of t h e m a t h e m a t i c a l

the flood wave propagation

has been

m o d e l of t h e A r n o in natural

has been built in an interactive

or

mode in order

to facilitate the setting up of the mathematical model and the simulation of the p h e n o m e n o n under different conditions. In

particular

w e shall try to emphasize h o w the use of an interactive

system and of a video unit has m a d e it very easy to use the model for any purpose. T h e flood routing problem lies in studying the propagation of a flood w a v e along the river, determining the variation of the wave's shape and height in its motion towards the sea. This is achieved by determining the evolution of the k n o w n water levels and diseharges at the initialtime (InitialConditions), according to the variation in time of the water levels and/or discharges at the starting section, and the inflow due to the tributaries (Boundary Conditions). The solution of this problem could then allow us to foresee the water levels and discharges in any section of the river at any time, assuming that the river bed geometry and its resistance to the water flow are known. The Mathematical M o d e l T h e solution of the flood routing problem in rivers needs the construction of a mathematical model, that is a set of elements and relationships, mathematical or logical, between them. In the ease of flood routing the mathematical relationships are the well known e q u a t i o n s (fig. l ) e s t a b l i s h e d which express these equations

mass

by Barr$

and momentum

[¥evdjevieh]

d e S a i n t V e n a n t i n 1871

conservation

implies

respectively.

substantially

[ De Saint Venant ] , T h e v a l i d i t y of

that the components

of t h e

232

water veloeity and acceleration, i.e.

other than the longitudinal one, are negligible,

t h e m o t i o n of t h e w a t e r c a n b e c o n s i d e r e d These

equations together with the levels z and discharges

u n k n o w n s of t h e p r o b l e m , superficial

contain geometrical

width B, hydraulic

q is the lateral

quantities,

resistance

Q, w h i c h a r e t h e

(like wet area A,

r a d i u s R), a n d o t h e r q u a n t i t i e s

n a t u r e of t h e r i v e r b e d , l i k e t h e p a s s i v e Furthermore

as unidimensional.

depending on the

J of t h e b e d t o t h e w a t e r flow.

i n f l o w o r o u t f l o w b y u n i t of l e n g h t a n d g t h e a c c e l e r a -

t i o n of the gravity.

.

.

.

.

.

B

.

v = mean velocity Q=A.v R = A/P

horizontal datum aQ ~x

+ R a__!_z- q = o ~t

az + j + 1 {av v av~ ax ~-\~ + ax/=° with

J= n2 v2 R,~

Fig.

The means

description

of the geometrical

of several

cross

of ±he watercourse

and

In the case example, Each

distance

is defined

Many formulas

shape

the number

on the accuracy

of our model

the mean

section

sections,

I

we

of the river of which

want

of the river Arno, between

by several

two

bed

depends

can be obtained

by

on the regularity

to achieve. [ Gallati

consecutive

and others

sections

] , for

is about

400 meters.

points of its boundary.

are available in order to estimate

the passive

resistance;

in

233

our m o d e l w e have adopted the M a n n i n g formula containing the coefficient n which depends on the nature of the river bed, and which, obviously, cannot be m e a s u r e d directly, F u r t h e r m o r e in watercourses there are other factors which cause energy losses, like, for example, bridges, bends, sudden broadenings or narrowings and so on. All these factors have not been taken into account separately, but they have been included in the above mentioned B/fanning formula. T h e coefficient n then m u s t be deduced by choosing the value which gives us the best fit of the k n o w n levels or discharges with the c o m p u t e d ones

[Gallati and

others] . In doing so, however, besides the friction losses, all the other approximations m a d e in the geometrical description and mathematical formulation of the p h e n o m e n o n also affect the value of this parameter; then it, like all the other p a r a m e t e r s which influence the p h e n o m e n o n ,

completely looses its original

m e a n i n g and b e c o m e s simply a fitting parameter. In order to solve the Saint Venant equations an original finite difference implicit s c h e m e has been adopted [Greco and Panattoni ] . W e

did not use an

explicit s c h e m e because of its constraints on the size of the temporal step, which preliminary investigations s h o w e d to be too little (a few seconds) for our problems. T h e use of an implicit scheme, although a little m o r e complicated, allowed us to use temporal steps of one hour, or more, which are m o r e suitable for flood routing problems. H o w e v e r it is noteworthy that, if w e use the N e w t o n m e t h o d to solve the nonlinear algebraic system, resulting f r o m an implicit s c h e m e , the matrix of the coefficients is always very simple

[ G r e c o and Panattoni] and this

m a k e s the solution of the s y s t e m very easy. Using this scheme, then, w e are able to c o m p u t e f r o m the k n o w n state of the river, that is, f r o m the level and discharge values in any section, at a given moment

T, the s a m e quantities at next time T+~T.Starting then f r o m the knowledge

of the river at the initial time, w e can determine the levels and discharges in any section at any time. In those sections, for which recorded data are available, w e can then c o m p a r e the m e a s u r e d and c o m p u t e d hydrographs and get s o m e idea of the model's efficiency. With the help of these c o m p a r a i s o n s w e can, as w e said before, set up the m o d e l varyingthe p a r a m e t e r s in order to obtain the best fit between the two hydrographs. A n d it is specially in this connection that the use of an interactive computing s y s t e m and of a video unit is very helpful. In fact the c o n t e m p o r a n e o u s influence of the various p a r a m e t e r s and their interdipendence are not k n o w n a

234

priori,

and by means

of a n i n t e r a c t i v e

meters

in an extremely

parameters

results

of h y p o t h e t i c a l

existing or in project,

The

Computation

System

m a k i n g t h e e v a l u a t i o n of t h e i n f l u e n c e of t h e

c a n b e o b t a i n e d a l s o d u r i n g t h e u s e of t h e m o d e l f o r t h e floods, for the management

of h y d r a u l i c

f o r t h e r e a l t i m e c o n t r o l of t h e r i v e r ,

the structure

of the Flood

things,

it assigns the default values to all the variables.

nicate,

by means

of a communications

terminal,

INITrAt

2. I

either

a n d s o on.

Routing Interactive

(FRISS).

The first step involves the initialization of the System,

Fig.

works,

System

The figs. 2. i and 2.2 show Simulation

of t h e p a r a -

on them very easy.

The same advantages simulation

we can vary the values

s i m p l e w a y a n d b y t h e u s e of a v i d e o u n i t w e c a n i m m e d i a -

tely display the complete various

system

and, among

The user,

with the system

then,

other can comu-

defining :

235

Fig.

2.2

~ PLOTTI~R

the watercourse

reach,

NOT

natural or artificial,

by means

of i t s s t a r t i n g

and

ending sections, - the geometry

of t h e r e a c h ,

either through the cross

sections which describe

it, or through its slope and shape, - the hydraulic

roughness

of t h e r e a c h (i. e. r o u g h n e s s

formula

and

coefficient), - the flood event, either past or hypothetical, - the boundary conditions, -

the l a t e r a l inflow or outflow,

- the numerical solution algorithm, -

the size of the time step,

- the accuracy of the iterative p r o c e s s , -

the use of the Jones correction formula for the stage-discharge relationship.

When the problem is completely defined, the program acquires the geometrical data from a river

geometry

file and then acquires

flood data from the past

236

flood events time

file or from

a video unit (hypothetical

2.2 shows

ges are computed At any time for instance -

-

the computation

a real

loop,

in which

the water

T il is possible

to display on the video unit some

profile along the reach,

the shape

of any chosen

Furthermore

cross-section

hydrograph

red with the measured

and the water

at any chosen

it is possible

and to restart the computation

gauging

to interrupt

from

select and change any wanted

unit, the plotter,

and the communications

3 which

profile of a river from

the printer

relating to the river Arno, has been

reach

the mouth.

obtained

about The

from

eighteen

lower

IBM

of the reach

-

to choose

compa-

are shown

FLOOD WAV~ 0~ HOUR

features

terminal. in the Figs.

long starting from

P~SA SC(~-NTtFIC CEN~ER W~,~'S

after storing

the outputs: using the

the video unit, represents

kilometers

~

and,

the processing

line is the river bottom

~qDPAr.,ATII~'IOF ~

3

section

T,

time.

it is possible

results,

level at that time

the computation,

At the end of the computation

Fig.

as,

one.

the situa±ion at the interr, lption time,

Fig.

pictures,

:

the water

Sovne

levels and dischar-

at all times.

- the computed

54 km

or directly from

system. Fig.

video

flood),

ARNO RI',,~'R

~ D~C i ~ 2.4

3 to 8. a longitudinal

the section

and the upper

one is the

237

tBH - P ~ A

SCIENTIFLC CENTER

FLOOD

WAVE O~

~B D~_C - ~ 2~

HOUR

Fig.

water

level as computed

by the model

Fig.

a specific

for a given time

of the flood of December

1948. 4 represents

with the computed Finally

water

in Fig.

of the considered

section,

level at a given time

5, besides

river

cross

reach,

about

40 km

for the same

the measured

flood wave

both the measured

F~D

Fig.

~

FLOCID I~VE~

WAVE ~F

5

"f"1~( ~ ' ~ )

HBUR

the mouth,

flood• in the upstream

section

(the dotted line) and computed

IBM - PISA SCI[NT~FIC C E ~ pR~Pt*,G~T~ ~

from

~R'~4B~VE~

~ DEC 1~4~

43

238

PIENA DEL

25 NOV 19~9

:1

L~

TEMPO-=,.

Fig.

( t h e f u l l line} w a v e s i n t h e d o w n s t r e a m f l o o d of D e c e m b e r

6

section are drawn.

1948 t h e r e c o n s t r u c t i o n

It i s c l e a r t h a t f o r t h e

of t h e a c t u a l p h e n o m e n o n i s q u i t e g o o d

and the differences between the computed and measured

v a l u e s of t h e w a t e r l e v e l s

never exceed a few centimeters. F i g . 6 on t h e o t h e r h a n d h a s b e e n o b t a i n e d f r o m t h e p l o t t e r unit. On it you can see, besides the upstream hydrographs

of t h e N o v e m b e r

measured

wave, some reconstructions

of t h e

1949 f l o o d a t s e v e r a l g a u g i n g s t a t i o n s b e t w e e n t h e

f i r s t o n e , s t i l l 54 k m f r o m t h e m o u t h , a n d t h e s e a . Finally in Fig.

7 and Fig.

8 similar

g r a p h s r e l e v a n t to t h e flood of J a n u a r y

1948 a n d t o t h e f l o o d of J a n u a r y 1949 a r e s h o w n . Conclusion T h e s e a r e o n l y a f e w e x a m p l e s f r o m a m o n g t h e g r e a t n u m b e r of e v e n t s w e n e e d e d t o t a k e i n t o a c c o u n t a n d to p r o c e s s

in o r d e r to a c h i e v e an a c c u r a t e k n o w l e d g e

of t h e r i v e r b e h a v i o u r d u r i n g f l o o d s ; a n d , i n c o n c l u s i o n , w e w a n t t o e m p h a s i z e a g a i n h o w t h e u s e of a n i n t e r a c t i v e s y s t e m a n d a v i d e o u n i t h a v e m a d e it p o s s i b l e t o s e t

~'~"OdN31

[ -_-j " ~ - ~

I .... ~LILL.................2SJTL 8 "I.~I

6h6! N30 EO

-

' ' ' 7 ' ' '

. . . .

'

. . . . .

'

. . . . .

'

. . . .

,

--

,

. . . .

,

. . . .

,

. . . .

,

. . . .

130 ~N31d

,

. . . .

,

"CY \17-,-s,~-7:

........

P, "II~ ~I

;

}

'

i

OhSI N3~ L~

7~0 ~Nl~a

6~"

240

up our model

very quickly,

through

the fitting of a grea% number

of past events,

and to use it very easily for any purpose,

such as studying the concurrence

different parameters

supplying the engineering

necessary engineering

on the phenomenon,

data (water levels and discharges), hydraulic

works

simulating

of

design with the

the behaviour

of

to be built, and giving objective and timely

data for

flood control.

References

De Saint Venant, B. "Th~orie du m o v e m e n t non-permanent des eaux avec application aux crues des rivi~res et ~ l'introduction des mar~es dans leur lit" Aead. set. Paris

Comptes rendus, V. 73, p. 148-154, 237-240

P a r i s 1871. Gallati, M., G r e c o , F . , Malone, U . , and P a n a t t o n i , L. "Modello m a t e m a t i c o p e r lo studio della p r o p a g a z i o n e delle onde di p i e n a n e t c o r s i d ' a c q u a n a t u r a l i " XIII Convegno di I d r a u l i c a e C o s t r u z i o n i I d r a u l i e h e , Milano, Sept. 1972. Greco,

F.,

and Panattoni, L. "An implicit method to solve Saint Venant equations". To be published.

Y e v d j e v i c h , V. "Bibliog:zaphy and discussion of Flood-Routing methods and unsteady flow in channels" Geological Survey Water-Supply, 1690. W a s h i n g t o n 1964.

Paper

AN AUTOFtATICRIVER ?LANNING 0PERATING SYSTE~ ,,,(ARPOS) Enrico Martins

-

Bruno Simeone

-

Tommaso Toffoll

IAC, Consiglio Nazionale delle Ricerche, Rome,Italy

Summar~ An analysis of the general structure of LP models of water resource systems and of the operations required in their constructio~ suggests the opportunity of designing a special purpose interactive system to assist the engineer in developing such models. As a working example of the ideas discussed, a simple operating system has been implemented. I. Introduction The optimal planning and management of complex water-resource systems requires the development of adequate models. Linear-programming (LP) models are almost invariably used as a starting point, even though they may be complemented,

at a later stage, by more

refined descriptions that use simulation or non-linear techniques (Buras and Herman, 1958, ~ a s s ARPOS is a special-purpose

and Hufschmldt,

1962).

interactive system (Galligani,1971)

aimed at assisting the analyst--typically a systems-oriented hydramllc engineer~-in the construction of ZP models of a water-resource system. In designing ARPOS, the development stage of a model has been kept in mind. In fact, models evolve together with the analyst's appreciation and understanding of the system's determining features. The model current version's inadequacies are a source of feedback onto earlier stages of t~ne construction process. Structural and numeric data are repeatedly reevaluated until a satisfactory behavior is achieved. On the other hand, since river models tend to be quite large (thousands of rows in the L P m a t r l x ) ,

every iteration of the revi-

sion process may require a large amount of error-prone work. For those reasons, aside from providing a general framework where

242

constructive activities and feedback loops are conveniently placed, ARPOS aimed at providing fast and reliable response to user feedback. This was achieved after an accurate analysis of the operations involved in the modeling process. In brief: (a) The whole process can be broken down into a sequence of welldelimited operations. (b) Many of the operations are well-defined and can be easily automatized. For instance, given a network representing the basin topology and the location of relevant activities (fig.l), the equations expressing continuity conditions at each node as linear comstraints can be generated automatically in a straightforward way. (c) Other operations, like, for instance, the removal of infeasibilitles, are required at certain well-defined stages of the process but are less amenable to a formal description and must rely ultimately on the analyst's judgement. Basin

network

~I

I R P A W E

~ R

~ .... i.....

P

= = = = = =

Inflow Reservoir Power Plant Agricoltural district Industrial area Ecological constraint

243

A more specific analysis will be carried out in Section 2. Section 4 summarizes the results of the synthesis process, in the form of an actual operating system, while section 3 illustrates some of the techniques on which such synthesis is based. 2. The ARPOS structmre. ARPOS takes into account the fact that certain features are almost invariably present in LP basin-planning models. These are usually multi-period,

which often causes the model to require thousands

of rows in the LP matrix. On the other hand, the model structure can be described more conveniently by considerably fewer symbolic constraints in which the time parameter appears only as a dummy index(+);" to this, a list must be added containing,

for all periods, only nume-

ric data. Such a time-independent symbolic-constraint

description

represents a convenient break-point in the modeling process: it has a familiar aspect for the hydraulic engineer, it simplifies the task of revising the model, and at the same time can be used as output or input of automatic procedures. Basic to the ARPOS structure is also the concept of standard constraints,

defined as those that represent continuity conditions

at individual nodes~ or flow bounds at individual nodes or branches. Standard constraints known models,

(a) make up a large part (90% or more) of most

(b) are maintained in successive versions of the same

model, and (c) can be automatically written down in symbolic form once the basin network is given. It is possible to express as standard constraints: - physical constraints:

continuity equations at each node;

- technological constraints: upper bounds due to the maximum capacity of channels, reservoirs,

and power plants;

i"+).... For instance, continuity at node 6 (fig.l) by the equation ~ + R6,t-R6,~+I = ~6,~÷ ~,t (I reservoir storage,'X outflow, P flow through the nes, all these variables beeing referred to node

can be expressed inflow, R current power plant turbi6 at the period t).

244

economic constraints:

-

lower bounds due to minimum local water

demand for agriculture, industry, or power generation; ecological constraints: lower bounds locally imposed on flows

-

for pollutant dilution, upper bounds for flood control; etc. The structure of the standard submodel (i.e., that consisting of standard constraints) makes it possible to test its feasibility by means of a special-purpose algorithm, namely, a variation of the well-known Hoffman-Fulkerson algorithm. Other constraints, which we may call peculiar (for instance, global economic constraints) vary from model to model and must be supplied directly by the analyst. Basically, ARPOS supports the following logic pattern of operations (fig.2): Figure 2

/

f ~

.)

t

L_

A c~vi~ie..S I

t

-

ffDaTacol~e~T~o~ ~ml~xeeh~alltor

5T~q~OI~RDS~J~OO~L

lI i

Basi~

el: 5UBI4ODEL)~---~ da'1'-a

1

I

I

c'r;~Tfai v',T$

I I \ \ 2.

245

It should be pointed out that ARPOS performs the feasibility test of the standard submodel at an earlier stage than the usual procedure

(Phase I of the Simplex Algorithm for the entire model).

This allows the user to enter the first ZP run with a higher degree of confidence on the model consistency. 3. Certaln techniques employed in ARPOS. 3 . 1 F e a s i b i l i t ~ test.

Checking feasibility of the standard subsi-

stem may be viewed as a problem of existence of feasible flows in the network E (+), which is defined as follows, starting from a given basin network B (flg.1): I) We produce p identical copies ~ ,..,Bp of B, where p is the number of periods, e.g., p = 12 for a monthly model. Thus to each node

x

of B corresponds a node

node x, we draw an arc from

x~

x~ to

in B~ . For each reservoir x~+ I (t = 1,...,p-I). This devi-

ce allows us to convert a network with memory into a memoryless network. 2) ~oreover, in order to get an all-activities-in-arcs representation, we add to each B~ a special node u~ from all nodes

which

and draw an arc to u~

are sinks or have some consumptive

activities; since our only concern is feasibility~we can always aggregate all the consumptive activities of each node into a single one. 3) It is also convenient, in order to have a simpler algorithm, to avoid sources and sinks. To this purpose, we add one more node v

and draw arcs from

each

u~

to

v

to all input nodes of BI, ..... ,Bp and from

v . It is convenient as well to eliminate multiple

arcs by inserting, when

necessary,

a fictitious node in some

arc. The network that is so obtained is, by definition, the extended network E (see flg.S).

,

(+)- Terminology is as in

~ord and Fulkerson, 1962.

246

,Extended network

for

FEAS_ --41-

"-1

1 I 1

-"

i

I I I 1 I I i I I

f

j,/

1

"

l

i

I t 1 l I l I I I t

J FIGURE 3

It is now straightforward to assign lower and upper capacities to the arcs of Eo These are defined: - for each reservoir arc (x~,x~÷ I ) to be zero and, resp., the maximum reservoir capacity - for each "consumptive" arc (y~,u t)

as the minimum, and resp.

the maximum requirement for the activity in node y in the period

t

- for any input arc (v,z E) both to be equal to the hydrological input in z in the period

t, and so on.

In order to check the existence of feasible flows in E, we use a known algorithm of Hoffman and Fulkerson (Ford and Fulkerson,1962 pag.52)o Without going into details, we only summarize the gross structure of the algorithm (called FEAS) in fig°4.

247

Macro

structure

o{

FEAS

A GUESS ~LO,,Z/ ) ,, -F = fo

A~~L:,gGE Y,IT :

~

,,,,Y

YfS'5.

CL~.N UP THe-

OLD ~AB~L~ 1 SBL~.CT 5OME A~t,~ INt=EA.~I BLE .__J

~ ~ E~kNT~L~OOGH ]

FIGURE 4 If the guess flow and the lower and upper capacities are all integral, it can be proved that, in a finite number of steps, the algorithm finds a feasible flow or detects the non-existence of such a flow. The assumption of integrality is not restrictive from a computational point of view. It should be

noted

that the algorithm is flexible in three

points: the initial guess flow, the selection of (s,t) and, inside the labelling procedure,

the selection of a suitable labeled

node

whose neighborhood is to be explored in order to create new labeled nodes.

We have implemented a version of the algorithm which ma-

kes use, in these three point~ of heuristics which rely on the special structure of E. So, for instance,

(s,t) is selected as down-

stream as possible in the hope that at each step many infeasibilities

will be reduced. The following advantages may be obtained With respect to the

normal procedure

(Phase I of the Simplex Algorithm):

248

a) a significant reduction in computation time, b) a better chance for the analyst to quickly discover causes of infeasibility, c) the possibility of obtaining, without making use of standard LP codes, policies which are feasible with respect to standard constraints, even if not optimal. 3.2 S~mbolic-constraint language and matrix generator.

This section

of ARPOS is described elsewhere (Toffoli, 1973) in greater detail. In brief: The input language to the LP-matrix generator should be as much as possible compact and user readable, especially since in ARPOS this language is also used as a vehicle for the user to add arbitrary additional constraints to the standard submodel. The hydraulic engineer's notation

for writing down water-system

constraints is immediately expressible in a context-free language defined by suitable productions. Moreover, such language is everywhere LRS (Look-Right-n, with n=1). In other words, in order to know what production generated a given syntactical element, it is sufficient to scan a_~tmost one symbol ahead in the input string. LRS-ness makes it easy to construct a fast single-pass parser for the input language. The formal definition of the symbolic-constraint language can be given in such a way that there is a strict correspondence between syntactic (Backus Normal Form productions) and semantic(matrix-generation operations) features. In this way, scanning of the symbolicconstraint list and matrix generation proceed jointly in an interpretive manner.

4. A s i m p l e O~erating syste m for ARPOS In order to provide a working example of the ideas discussed above, a Simple Operating System (SOS) for ARPOS was built in its entirety. In spite of its small size (it works well on an I B ~ 1 1 3 0 with 16K and one disk), SOS is fully able to support the construction of a model up to the point where standard ZP codes take over.

249

SOS has been extensively

used in developing LP models

for the

Tiber basin. SOS consists

of a set of ~rocedures

The stream of commands the interpreter programs

generated

into a sequence

is dynamically

(procedures)

selected

called by an interpreter.

by user is translated of procedure

constructed

calls.

on-line by

Thus,

an SOS

by joining together modules

from the following repertoire:

Es

Ex

LOAD - Load input deck in input data file. REVS - Edit or update input data file. LIST - List input data file. L~

- Lump numeric data for ~ . ~ elementary into data for m seasons.

time intervals

m

NETW - Generate PLOT - ~ p

river-system

river-system

STAT - Give statistic SYk~B - Generate

network.

tables

symbolic

network,

on network and activities.

equations

EDIT - Edit symbolic and objective

equations. function.

A L G B - List symbolic

equations.

FEAS - Test feasibility

for standard

Add additional

constraints. constraints

of current version of the

COXP - Compile LP matrix for input to ~£PS.

submodel.

250

A warning is given if a meaningless combination of modules is attempted. In this way, commands can be added or repeated with different data, in an interactive manner, and at any moment, on the basis of the information currently available about the model's status, the user will be able to revise previous operations or proceed with the construction.

References (I)

Buras, N., and Herman,T.:"A review of some applications of m~thematical programming in water resource engineerings"Prog.Rep. n.2, Technion, Haifa, Jul 1968.

(2)

For d jr, L.R.

and Fulkerson, D.R.:"Flows in networks~'Princ.

Univ. Press., Princeton 1962. (3)

Galligani, I.: Un sistema interattivo per la matematica humerica, Riv. di I nf0rmatica, vol.2, n.1, suppl. Apr 1971.

(4)

Maass, A., Hufschmidt, M.A.,and others:"Design of water-resource systems,"Harvard Univ. Press, Cambridge, 1962.

(5)

Toffoli, T.:'~ Precompiler for MPS,"IAC, 1973.

Acnowledffements The authors are indebted to Prof. I. Galligani for his suggestions and to Mr. A. Bonito for having written some of the programs.

ON THE OPTIMAL CONTROL ON AN INFINITE PLANNING HORIZON OF CONSUMPTION. POLLUTION, POPULATION AND NATURAL RESOURCE USE AlainHaurie Ecole des Hautes Etudes Conmmrciales and Ecole Polytechnique Montreal, Quebec, Canada

Michael Po Polis Ecole Polytechnique Montreal, Quebec, Canada

Pierre Yans ouni Ecole Polyte chnique Montreal, Quebec, Canada

ABSTRACT In this paper a five state variable economic planning model is presented. A Malthusian hypothesis is employed which gives rise to a zero-growth argument. The asymptotic behavior of the model is studied.

Classical results involving Turn-

pike theory are used in conjunction with recently published results on the infinite time optimal control problem to show the convergence towards a Von Neumann point of the optimal trajectory on an infinite planning horizon° 1. INTRODUCTION Recently J.W. Forrester in "World Dynamics" [l] proposed a simulation model to study the interactions between population growth, the use of natural resources, pollution and consumption°

The simulation results showed the catastrophic

effects which would result from a policy of "laissez-faire". been widely cormmented upon and criticized [2], [3].

ForresterTs work has

Two of the most repeated cri-

ticisms concern the high level of aggregation of the model ("five state variables to describe the world?"), and the Malthusian philosophy which gives rise to the fundamental hypotheses upon which the model is based°

Although these are valid cri-

ticisms they do not negate the utility of the model even though the conclusions drawn from the simulations may be questioned. Independently of Forrester, economists have developed a theory of optimal economic growth based largely on techniques of optimal control [4] -[7],

Effecti-

vely the complex phenomena of interaction simulated by Forrester can be studied under a form similar to that used by the economists°

Further, optimization techniques

may be used on the aggregated models to establish a direct relation between the conditions of the problem and the form of the solution which must be sought. In this paper an economic planning model having the same level of aggregation as Forresterts model, and under the Malthusian h~-pothesis is presented.

The

techniques of optimal control on an infinite time planning horizon are used in conjunction with the model to provide information on the nature of the solution which rmast be found°

Thus, a five state variable model of economic growth is proposed to

study, at a macro-economic level, planning policies concerning consumption, the use of replenishable and non~replenishable natural resources, pollution enmuission and

252

population growth.

Such a model can be useful in two ways:

(1) A price mechanism sustaining the optimal

policy can be defined and eventually

an optimal fiscal policy leading to optimal or suboptimal paths can be proposed (See Arrow and Kurz [8]). (2) The asymptotic behavior of the optimally controlled economy can be studied. In this work only the second point is considered. gives rise to a zero-growth argument.

A Malthusian model

In the context of an economy without marked

technological progress the objective of zero-growth has already been proposed [9]. The planning model studied in the following paragraphs permits the rational choice between various equilibrium states which are all characteristic of zero growth° Further, these equilibrium states indicate how the capital structure of the economy can be transformed in order to attain a desired objective following an optimal trajectory consistent with the objective° The principal theoretical result presented here concerns the convergence to a steady-state-zero-growth-economy (¥on Neumann path) of an optimal trajectory on an infinite planning horizon°

This result is established by using a definition

of optimality recently proposed by Halkin [10] as well as classical results concerning the well known "Turnpike" theorem for a multisector economy [II], [12]. 2. THE SYSTEMIS MODEL Models used in the neo-classical theory of economic growth are highly aggregated.

Most often only one state variable is considered: the capital stock;

labor is an exogeneous variable assumed to grow exponentially with time~ and production is generated by the combination of labor and capital stock. In order to account for such mechanisms as pollution, population control~ extraction and exhaustion of natural resources, a model is considered with the following five state variables: K

A

capital stock.

Xl

=a non-replenishable resource i.e. minerals~ oil, etc.°.

X2

A

P

=a pollution stock°

L

~

replenishable resource i.eo forests, fish stock~ etCooo

total labor force.

It is also assumed that direct manipulation of the following economic variables~ is possible : Y

=

production level in the economy°

C

~

total consumption°

253

N

~

part of the revenue allocated to control the birth-rate.

R1

extraction rate of the non-replenishable resource°

R2

~ & =

extraction rate of the replenishable resource.

Q

~

level of economic activity dedicated to the regeneration of X 2 .

G

~

pollution emission rate.

These in fact are the control variables in the model°

The relationship between

state and control variables is imbedded in the state equations describing the dynamic behavior of the economy°

-

f(

=

Y -

~l

=

- RI

C

- N -(~K

(2.1)

(2.2)

X2 = 1~ (X2' P' Q) - R2

(2.3)

~'

=

0-vP

(2°4)

L

=

D(L, N)

(2.5)

Equation (2.1) represents the capital accumulation in the economy;

G

is the

capital depreciation rate. -

-

Equation (2.2) represents the extraction of the non replenishable resource° Equation (2.3) represents the extraction of the replenishable resource;

~(o)

is the rate of regeneration of X 2 o -

Equation (2.4) represents the accumulation of pollution.

v

is a natural elim~-

nation rate of pollution by the environment° -

Equation (2.5) represents the population dynamics°

D(o, • )

is the population

"growth" (it can be negative) rate. To a certain extent the mathematical description of the mechanisms involved in the economicsectors modeled by equation 2.1 - 2o 5, is inspired by the work in

Ref°

[13] - [20].

A more detailed analysis of the mechanisms involved in the manipulation of the control variables is necessary

to underline the physical constraints res-

training the set of possible choices for the control variableso

Stating again the

neo-classical hypothesis that economic activity is generated by the combination of capital and labor allocated to each sector; it is assumed that a fraction of the capital

K

and of the labor

L

is allocated to each of the areas of economic ac-

tlvity defined by the state equations.

Each of these fractions is identified by a

subscript referring to the respective control variable.

254 They are : ~

KI, K2, KQ, KG

~, h, h' LQ, % Two physical constraints arise inmmdiately, they are: ~+

~+

+

h

(2.6)

K 2 + KQ + K G • K

+ L2 + LQ+ L~< L

(2.7)

Denote the maximum~ obtainable rate of "activity" in tible with the allocated fraction of ~.

K

and

each sector,

L ~ by the variables

compa-

Y~ RI~ R2~ Q,

A set of "technological" constraints arises from the fact that the preceding

variables constitute upper ~ bounds to the control variables. O
Io

f ( x 2, u 2)

with

u2 ~0

[O,l] , ~ u ~ > u

, %(x2, u2)~O

, ~j(~xl + (l - 4) x 2, J ) > o

such that: ~ll + (i - ~) 3 2

=

f(~x1 + (1 - ~) x 2, u ~)

This result is e a s ~ y demonstrated with the assumptions that

@(o)

and

D(.)

are

concave, and further, that they are linear with respect to the state variables. Without the linearity assumption no demonstration has yet been found. general functions ~(o)

and

D(o)

For the more

the result stated above has to be considered a

basic assumption to which can be applied the comments contained in Remark 2.1. The economic model is to be completed with a welfare criterion defined on an infinite planning horizon.

It is assumed that the objective of the economic po-

licy is to maximize a functional: J where

W(o )

choice of

ki > / O

f

W(C, L, P, G, XI, X 2) dt

is restricted to be a concave function of its argun~nts.

W(o)

is interesting because of its simple implementation:

J with

=

--

and ~ X i~< i o i

fo° (XlC + k2X1 -

X3P - X4L) d t

A particular

256

Assuming convergence to the Von Neumann path, (see section 3) the equilibrium solution of the optimal control problem is parametrized as a function of the weights

ki

attributed to the respective variables in the welfare criterion.

Sol-

ving the algebraic equations defined simultaneously by the necessar~~ conditions for optimality (on the infinite horizon, Ref. [ I03) and the equilibrium (steady state) condition, is a relatively simple task. t The solution yields in fact the equilibrittm levels of consumption, population, pollution~ and the remaining stock of non replenishable resource.

This procedure gives the n~ans of evaluating the long term

consequences of a given choice of values for the weights

ki o

Furthermore a sys-

tematic exploration of the domain of equilibrium solutions as a function of the pal rameters X i will give a basis for assigning values to Xis compatible with the C long range objective of the econon~-~ i.eo expected per capita consumption ~ , population and pollution levels allowed, and amount of resource

X 1 to be protected.

Remark 2o3: The use of a weighted performance criterion may also correspond to the search for Pareto-optimal solutions arising from a vector valued criterion° deration of a non-scalar

Consi-

criterion is indeed most likely to occur in an actual

planning situation° 3o OPTIMALITY ON AN INFINITE TIME HORIZON AND CONVERGENCE TOWARD A STEADY-STATE Consider the d)~aam_ic system:

f(x, u)

=

x(O)

=

x°

(3.1)

given

x C X c E n, u¢ U ( x ) c E m f

is

C1

utility

with respect to xo

x

and continuous with respect to

=

0

where

w

is

u o The accumulated

satisfies :

C1

with respect to

x

w(x, u)

(3.2)

and continuous with respect to

u .

Definition 3.1: The control

~" : [0, ~] * E m

is optimal if the following conditions are

s atis fied : (i)

it generates a trajectory

( ~ , ~o ) : [0,

~)

+ 1 ~ En

such that

This is equivalent to solving a standard mathematical programming problem.

257

: EO, ~) * ETM

(ii) V ~

generating Vtg(t)

then

such that

(~, "~o ) : [ 0 , ~) ~ En + 1

~ U[g(t)] , ~(t)~ X

V¢>oVt >0 ~ T>t

such that

~o(T)< ~ : ( t ) + ¢

(3.3)

This definition of optimality on an infinite time horizon has been proposed by Halkin [10]o The following assumptions are inspired by the Turnpike theory [II], [12] in mathematical economics.

Assump,tion 3.1: The set ~ ~ {(x, ~, ~o ) : ~xcX , ~ ucW(x) =

f(x,

u) , ~o

= w(x, u ) }

~

such that:

o l o s e d and oonvex.

Assumpt,~Rn 3.2: If

(x, ~, ~o)eQ

and

x' > x

, ~ _ ~ Xo ~ x'~< ~

then

This assumption is equivalent to the free disposition of goods. Assumption 3.3: The following propositions hold

-~x~X,3u~U(x) -V~>o,

-~c

sot.

s.t.

3 ~(~) s.t° ILxll>c

f(~, u ) > O

(3.4)

llxll< ~ ~ Vueu(~), Ill(x, u) , w(x, u)}l< ~(~) (3.s) d I/xl12 = < *Vu~(~),-}~-

The first main result can thus

x, f(x, u ) > < O

(3.6)

be established°

Lenmm 3.1: There exists

(i)

4~ w cR

and

p"-¢En ~ p'~ ~ 0

such that

(3.7)

w* =fMax Xo : (x, ~, ~o)¢Q , ~ > 0 7

(3.s)

II Proof:

Define the sets:

V ~

{(~, ~0 ) : ]xeX

s.t.

(x, £, ~o)¢Q]

v+ ~ f~, ~o)~v : ~ 1>o]

For this assumption to hold Equation (2.4) must be written

(3.9)

(3.1o)

- P = - G + vP o

258 It is clear from Assumption 3.1 and Equation (3.4) that convex~ and that from Assumption 3.3, Equation (3.6)

w

(x, Xo)OV÷l

o

can be d e f i n e d on If

V+

(~q"-, w~) cV÷ , t h e n Since

V

is non-empty closed is compact.

Thus:

:

which establishes (i)o

not, it would be possible t o find impossible.

=

V

V+

(~'-, w~"-)¢~V , t h e boundary o f the s e t ~0

and

w>~*

is closed convex there exists

s.t.

(](, w) ~V +

(W, p) e En + i

p>0

If

which is

such that:

V(~, ~o ) evnJ"- + < p , ~e> >n~o + < P ' ~ > FromAssumption 3o2~

V .

(3.n)

and it is always possible to choose x = 0 o Thus % p is the corresponding value of p .

is positive and can be choosen equal to one;

Equation (3.8) is then a direct consequence of Equation (3.11). Definition 3o2: Following the economistTs terminology a "point of Yon Neumarm" is a vector

(x, ~, % ) ¢ 0

such that ~o + < P

-k ~' ~ > =

w-,,~

(3.12)

The system defined in (3ol) and (3.2) is said to be regular if the set of ¥on Neumann points has only one element:

(d ~, o, ~ ) Remark 3.13 w state,

is also the maximal utility flow which can be maintained at steady(x~, ,2'-) is the solution to the following optimization problem:

i.e,

( Max w(x, u)

I

f(x, u) = 0

ucU(x) x c X

Now c o n s i d e r t h e convergence towards a Yon Neumann p o i n t o f an o p t i m a l t r a j e c t o r y on an i n f i n i t e

time horizon°

Once a g a i n a c l a s s i c r e s u l t

from t h e Turnpike t h e o r y

[ 11] i s usedo Lemna 3.2: Let

F

y =~ (x, ~, ~o)¢

be the set of all Von Neumann points, define t h e distance: d ( y , F)

~

M i n l l y - zll $¢F

F C E 2n + 1, for

259 Thus : V¢ ~>0, ]8(e)

Ilxll
t l ~ ( t )

= u* s.t.

f(~"-, u*)

V t > / t I xo(t)

= xo(t 1) + ~ ( t

=

o

Thus : ~T

- t l)

(3.15)

260

Now c o n s i d e r t h e o p t i m a l t r a j e c t o r y ;

~ " : [0, ~] * En o

Since

II x°I1 < C

and from

assumption 3°3, E q u a t i o n ( 3 . 6 ) : ..9 llx(t)11O

Thus, from Lemma 3°2 and Equation (3.8) the following must hold:

~j(t)

2

=

w [ ~ ( t ) , ff~(t)] dt

~< tw* - < p * ,

( ~ * ( t ) - XO) > - 8 %

(3.16)

whe re :

a ~[{Te [O,t] : d(y*(~), F ) > ~}] Pt =

(3.17)

I f 3°14 does not hold, then tl~im pt = m , and thus (3~15) and (3.16), denoting the largest component of

p*

as P~ax ' imply: ] V ¢ > O ~T

~:(t')

s.t. V t ' >

+ ¢.n'

iff

SOME THEOREMS. I. Blackwell Theorem. n is more informative than ~' iff n is quasigarbled into n'. Corollaries:~If n' is coarser than n, then n is more informative than n' ~ T h e lattice induced on the set of information systems by the relation "more information than" is identical with the lattice induced by the relation "quasi-garbled into". 2. Some "desirable" properties of concave symmetric function§ on the probability space, un(ql,...,q n) has a unique maximum on the n-simplex, and by symmetry 2.1

un(~,., ",~ i) ~

un(ql , .

. .

,qn )

;

268 and since [~+1 (1

I

~,...,~,o1

:

un(1

1

K,...,~)~

it follows that 2 2

U n+l

•

1 --l.1 )> (n-$-~'''''n+l

un(1 ... 1

un(~,, "''ni)

1 1 U m (~,...

;

and

if n>m.

(Note: 2.1 and 2.2 hold also for quasi-concave symmetric functions).

8. Non-negativity of "uncertainty remoyed": for all Ug{U},

u(z) - u(zIY)

~

0,

("added evidence cannot increase uncertainty").(Note:

3. holds for

all concave functions. See De Groot.) 4. Comparative "informativeness" vs. comparative "information transmitted": n>n' iff U(Z) - u(ZIY) ~/ u(z) - U(ZIY') (Note: 4. holds for all concave functions.

for all UE{U}.

See Marschak and

Miyasawa .) Corollary: Suppose H(Z) - H(Z[Y)> H(Z) - H(ZIY'), but there exists U~{U} such that u(z)

- u(zI~)

< u(z)

- U(zIY').

Then there exist ~,~ such that V ~(n) < V ~(n').

5. Properties exclusiv e to e ntropy. H. 5.1

Additivity. U(plql,...,plqn,p2q2,...,pmq n) = U(p) + U(q)

iff U = H.

(This permits to "measure information". What for?)

269

5.2 Minimum expected length of d ecodable qi = probability

code word.

(prior or posterior)

of the word, encoding expositionjneglect

Let

of event z i and I i = length

it in alphabet of size r. For the sake of

the integer constraint

on it, a constraint

which makes the result an asymptotic one (Shannon; Feinstein,

Wolfowitz).

should be decodable.

Ii r is sufficient

But retain the constraint

see also

that the code

For this, the "Kraft inequality"

-i i

~

1

and necessary.

Write 1 = vector[it] ,

--l.

L = {i: [r

lE

1}. Then

min I&L

U =

[qili = U(q)

iff

[qi l°gr (I/qi) ~ H

since, by simple calculus,

the optimal

code word length

I i = logr(i/qi). The result agrees with the characterization measure of disorder

if high disorder

a large number of parameters rent appears

to describe

to me the mathematical

a state.

derivation

the logarithm of the Stirling approximation of a given frequency

COMMUNICATION

distribution:

of entropy as a

is understood

to require Quite diffe-

of entropy via

of the probability

see, e.g.,

Schroedinger.

COST, OTHER SYSTEM COSTS,

AND SYSTEM VALUE.

The cost of communication late long sequences

-- of s±oring the messages

encodable

economically,

transmit and decode them --, does depend cally) on the entropy parameters distributions.

to accumu-

then to encode,

(at least asymptoti-

of the pertinent

probability

But these steps of the communication

process

form, in general,

only a part of the information

already described

such a system as a chain or, more generally,

a network of (stochastic) symbol-outputs.

transformers

The communication

system.

of symbol-inputs

We have into

steps are such transformers.

270

But they are preceded by, for example, (generally unobservable, into observed

storing,

of data into statistical

possibly to be communicated encoding,

transmitting,

estimates

or pre-

to the decision-maker

decoding.

cost is in general dependent

via

Each transformer

associated with the cost of acquiring and operating operation

of

events

data, related to the events by a likelihood matrix;

and the transformation dictions,

the transformation

e.g., future) benefit-relevant

is

it, and the

on the transformer's

inputs

and is therefore random. There is no reason why thecosts of sampling or of manipulating data ~ for example~ should be related to any entropy-like matrices

expression

involving the relevant

and the input-probabilities.

information

The expected

likelihood cost of The

system q will depend on the prior distribution

~ and

will be denoted by K (~). If we neglect as before the costs of decision-making,

the problem of the chooser and user of the system

n is to maximize

the expected net benefit,

between the expectations max n

max

where Vws(n),

of benefit and of cost:

(B 5(e,n)

= max n ~Vw~(n)

i.e.~the difference

- K (n~

- K (n))

as defined before,

,

is the value of the information

system ~. We have shown that this value~ too, like the system costs other than those of communication, parameters.

If the assumption

is dropped, max

the problem becomes: ,n (B 5(s~n)

where C (e,n), the expected probability

does not depend on any entropy

of equal costs of all pure strategies

- Kw(n) - C (~,n~

,

cost of strategy e~ depends on the

of its input y, and hence on the given ~ and the chosen

n. Now it becomes

impossible

to separate the expected benefit

B ~(s,n) when maximizing with respect to e. The "value of the information

system"~

mixed strategies, e = [ ~ y -~ ~.

V ~(n) loses its significance. expressed

by,generally

where eya=P(aly) , must now be considered.

is no reason to relate strategy

But again, there

costs to any entropy formula

I feel therefore that the assertions justified.

At the same time,

"noisy", Markov matrices

made in the SUMMARY are

271

BIBLIOGRAPHY Acz~l, J. on Different Characterizations of Entropies. Probability and Information Theory, 1-11. Behara, M. et al., eds., Springer (1969) Blackwell, D. Equivalent Comparisons of Experiments. Ann. Math. Stat. 24, 265-72 (1953) and

Girshick, A. Theory of Games and Statistical

Decisions, Me Graw Hill (197o) De Groot, M.H. Uncertainty~ Information and Sequential Experiments, Ann. Math. Stat. 33, 4o4-419 (1962) Optimal Statistical Decisions. Me Graw Hill (197o) Feinstein, A. Foundations of Information Theory. Mc Grsw Hill (1958) Marschak, J. Economies of Information Systems. J. Amer. Star. Ass, 66, 192-219 (1971) Optimal Systems for Information and Deoision, Techniques of Optimization. Academic Press (1972) and Miyasawa, K. Economic Comparability of Information Systems. Intern. Econ. Rev., 9,137-74 (1968) Savage, L.J. The Foundations of Statistics, Wiley (1954) Schroedinger, E. Statistical Thermodynamics~ Cambridge Univ. Press (1948) Shannon, C. The Mathematical Theory of Communication. Bell Syst. Tech. J. (1948) Wolfowitz, T. Coding Theorems of Information Thgory. Springer (1961)

ON A DUAL CONTROL APPROACH TO THE PRICING POLICIES OF A TRADING SPECIALIST

Masanao Aoki Depar~mmnt of System Science University of California Los Angeles

I.

INTRODUCTION

We consider a market in which there is a specific economic agent or authority who sets the price, and trading takes place out of equilibrium. keteer following Clower [8].

We call him a mar-

A trading specialist is one example of a marketeer

whose pricing policies are the main concern of this papem.

We assume further that

he does not know the exact demand and supply condition that he faces in the market, i.e., he has only imperfect knowledge of the market response to a price he sets.

He

does have some subjective estimate of the market response as a function of the price he sets. See Arrow [2] for a related topic. Denote by f(p;8) his subjective estimate of the market response (for example, excess demand for the c o ~ t y ) ,

where 8 is an element of a known set e . In other

~rds, in the opinion of the marketeer {f(p;@), 8 ~ @} represents a family of possible responses to his setting p . By specifying 8, a specific response is chosen as his estimate.

Assume that the true response (unknown to him) corresponds to f(p;8*)

+ ~ , 8* g 8 ~ where ~ is a random vamiable to be fully specified later.

We take 0

to be time-invariant. Therefore~ the actual response may deviate from f(p;8) by: (i) the agent's estimate of the parameter e, being different from the true parameter 8*; and by (ii) the random variable ~ which may be used to represent the effects being Lm]q%own to be systematic to the agent.

Various anticipated or systematic trends % such as price expec-

tation on the part of the buyers could be modeled. The probability distribution of is assumed to be known to the agent.

The type of consideration to b e presented below

can be easily extended to the case where the distribution is known imperfectly; for example, up to certain paran~tem values specifying the distributions uniquely.

This

added generality is not included in the paper, since it represents a straightfor~4ard extension of this paper. %To consider these~ p in f(p;e) must be replaced by the history of past prices. do not consider this case explicitly in this paper.

We

273

The model to be discussed in this paper could arise in several economically intere&ting contexts ~ for exan~!e ~ in considering optinml price setting policy of a /~onQpQlist ~, (Barro [5]) or a trading specialist or inventory-holding middleman in a flow-stock exchange economy.

For a non-random treatment of a pricing question

under imperfect knowledge, see Hadar-Hillinger [10]. In this paper we present the case of a trading specialist by assuming that the agent holds sufficient inventory so that he can complete all transactions.

The non-

price rationing behavior will not be discussed. In Section 2, we formulate the model.

The process by ~ c h

the marketeer's sub-

jective estimate is updated with additional observation is discussed in Section 3, and the derivation of optimal and suboptimal (second-best) pricing policies are derived in Section 4.

The last Section discusses some extensions and approximations.

2.

MODEL

Consider a trading specialist dealing in an isolated market who trades in a single good.

As pointed out in [7], monopolistic or oligopolistic price adjustment

can be formulated in an entirely analogous manner.

The only main technical differ-

ence is due to the fact that the price and output rate must be treated as decision variables, rather than just the price. We consider a Marshallian "short period" market.

At the beginning of each tra-

ding pemiod (day, say), he posts the price at which the good is to be traded. are several possibilities.

There

His selling and buying price will be different by a fixed

percentage, or he faces i~recisely known demand but known supply, and so on.

In

each case, it may be reduced to a situation with his subjective estimate of excess demand being given by f(p;8) with appropriately chosen e . Since he trades in disequilibrium, he is not certain what current price p will clear the market.

When he

changes the price of the good, he has only his subjective estimate of the effects of the price change.

Faced with this uncertainty, however, the trading specialist sets

the price which, in his estimate, will maximize some chosen intertemporel criterion function. We assume, therefore, that the trading specialist's criterion function over the next T periods is a function of xt, t = 0, ..., T where ~

is the excess demand at

time t xt = d t - s with ~

t

and st being the market demand and supply at the t-th marketing day. Note

that x t becomes known only at the end of the t-th market day.

Barro's approach ~ay be considered to be a special case of this paper where @ is a singleton. We do not consider price adjustment cost, however, in the criterion func~ tion of the specialist. This is done for the sake of s~nplicity of presentation.

274

We examine an optimal price setting policy by formulating this problem as a parameter adaptive ( ~ a parameter learning) optimal control problem, since he learns about the ~ o w n

parameters which specify the excess demand function.

We then apply

the dual control theory to derive the equation for the optimal price policy [3]. This paper therefore represents one possible mathematical formalization of price setting mechanisms discussed by Clowem [6-83, who emphasized the subjective nature of economic agents' decision making processes. Let St be the stock level of the conmDdity at time t . We have

%+1 : st - xt

(i)

Rather than assuming that the marketeer chooses p by setting f(P;0t ) : 0 , where 0t is his c ~ t

estimate of 8, we consider optimal and suboptimal pricing

policies by minimizing a certain criterion function explicitly. The criterion function to be minimized is taken to be T

J : k(ar+l) + [

t:0

where

Ct

~t(xt)

(2)

is a function of x t , and represents the cost of being out of equilibrium

(non-clearing market days). the end of the o ~ t

The function k evaluates the cost of the stock ST+1

at

planning horizon deviating from a desirable stock level. One

possibility is to set

k%+l ) : (ST+l_ ST l)2 , with __-~+Iplanned d e s ~ these two costs.

:

C2')

stock level; I is the parameter giving relative weights of

The time discounting factor can easily be incorporated but set

equal to zero here for simplicity. Another possible criterion function at the s-th market day is T

f-

T

to express the fact that the specialist has a preferred level of stocks throughout the planned period, and to deviate from it in either direction is costly. These criterion functions incorporate explicitly costs due to non-clearing markets and to deviation of stock levels from some desired levels, reflecting such considerations as perishable goods and storage costs. A third possibility for J is to include cost due to price changes.

The cost in

this case is due partly to the price stability and predictability which is relevant in minimizing search cost~ and not so much due to the set-up cost of price changes [13. This additional cost is not incorporated here in order to focus our attention on the parameter learning aspect. paper to be specific.

The criterion function (2) will be assumed in this

275

The criterion functions (2) express the possibility that the marketeer may wish to increase or decrease the terminal stock level by purposely setting price that makes his estimated excess demand flow at non-zero values. The stock of the c ~ i t y

held by the marketeer at time t, St, is therefore

given by t St+l : So - ~ XT(PT) T=0

'

where the excess demand x t is assumed to be a function of Pt only. We assume the market excess demand function xt(p) is to be given paremetrically as xt(p) = f(p;8*) + ~t

' 8. E 8

(4)

where 8* is the parameter vector that specifies a particular excess demand curve out of a family of such curves~ and where { ~ } is taken to be a sequence of independent~ identically distributed random variables %

with

E~t = 0 , a l l t >_0 var ~t = 2

< = ~ all t > 0 .

In this paper we assume ~ to be known. One example discussed later considers a linear excess demand function xt(Pt;~'8) = "~Pt + 8 + ~t ' ~ > 0 , 8 > O .

(4')

The unknown parameters are e and 8 %f in this case. It is assu~ed that the trading specialist has a priori probability density function for 8, Po(8) . The problem may now be stated as: Determine the sequence of prices Po'Pl ~ "" "~ ~ E[JIH t]

such that

t = 0,i, ..., T

are minimized subject to the dynamic equation St+1 = St -xtCPt)

t = 0,i ..... T

where

%Serially cor~lated ~'s can be handled with no conceptual or technical difficulties. The independence assumption is chosen for the sake of simplicity. %#Theoretically, there is no difficulty to increase the number of unknown parameters more to include c , say. See See. III. 2 of [3]. The computation becomes more involved, of course.

276

"t = {Xt-lPt-l' ~t-1 l %

= { % and a priori knowledge of 8}

As stated above ~ we do not impose the constraint St > 0 for all t explicitly, assuming that % is sufficiently large. At time t, St is a known number.

However, St+l~ .... % + 1 are random variables.

3. DERIVATION OF CONDITIONAL PROBABILITY DENSITY FUNCTION OF At the beginning of the t-th market day, the specialist knows the past and current stock levels, ST, T = 0,1, ..., t and from these past excess demands xT(PT) , T = 0 ,i, ..., t-i ; and he computes the posterior probability density function of 8, p(elHt) by the Bayes rule recursively from po(e) :

P{e IHt+z)

=

p¢elHt) p{:~t I Ht,e,p t) P(Xt IHt ,Pt )

=

p(e I~t) p(xt le,pt)

(6)

P(Xt IHt ,Pt )

where

p(elHo) :

p(xote,po) po(e) f p(xole,p o) po(e) dO e

where we compute p(xtle,p t) from our knowledge of the probability density function for the noise ~ . The trading specialist's uncertainty or his knowledge about the unknown market excess demand function at time t is sursremized by p(elHt). As a function of Ht, he sets Pt for the market period t . For example, if {t

is Gaussian, with mean 0 and standard deviation q, then we

have

PC~tle'pt) " 2/~i a exp_ 2~1 2 [xt_ f(e,pt)] 2 Unless po(e) is such that it admits sufficient statistics, the size of Ht linearly grows with t . The calculation involved in updating the posterior density function by (6) may be very large.

Only a few density functions have the pleasant propezi-y of '~epro-

ducing" the functional form after the transformation (G). The normal density function and Beta density functions are two well-known ones [ 3]. To illustrate, suppose that f(p;e) is linear as given by (4'), with e = Suppose further that po(e), is given by the normal density function with mean vector 8o and the covariance matrix Ao . By the repeated application of (6), we obtain the normal posterior density function p(81Ht) with:

277 ~tctxt+ ~tAo~o

(7)

where

xt =

=

[i .... i]

,

Ct

t-i 2 X1 + ~T=0 PT

'

t-i X2 + ~~=0 P

t-i 2 + ~T=0 PT

'

13 + t

where

I

and where X1 % X3 are defined by °2A-lo =

J

I X1 X2] t2

X3

It is instructive to write (7) as \

~t = ~I - Kt) t-1 ÷

A

t-i~[pt-i~ i ! xt-i(Pt-i

(7')

where

7bt-~ ~÷

cpt

-

~,:~)At-~/~t ~! ,;2kl

which shows explicitly how the latest piece of data xt_ I ~t-1 ) influences the trading specialist's estimate of 8 at time t, 8t " A similar, slightly more involved expression obtains if the Beta density function is used.

This function may be suitable to be used for ~ in the linear excess

demand function, since it is quite likely that 0 < ~'< i .

278

4. DERIVATION 0£ THE OPTIMAL PRICE POLICY In this section, we derive a sequence of optimal prices. ments in Sec. III of Ref.[3].

We follow the develop-

Last Sta~e Suppose Po'Pl' "" "' PT-I has been determined, and only chosen. PT is chosen to maximize

The excess demand at time T, x T , hence

ST+1

PT r e d s

to be

is predicted by

~'(kCST+l)I~T) : I k(,%+p P(ST+lIST,e 'PT) p~8 IHT) dST+ld8 (8) = I k ( ~ - f(PT;S)-() p(~) p(elHT) ded( where

p(~)

is the known probability density function for

~ , and

E(XT(PT)IHT}= I f(PT;e)P(elHT)de When k(.) and f(.) are specified, (8) can be evaluated explicitly, at least numerically. The optimum PT

is determined by

Max Elk(ST+I)+ ¢T(XT)IHT] . Denote the maximum value by Next to Last

YN "

Stage

The prices

PT-I

end PT are to be chosen to maximize

The computation of E[¢T_I(pT_I) IHT_I] can be carried out similarly

tO the

last

stage case. The other two terms are computed as

At time T-l, ~-i' ~-i' and ~ - i are known. XT_i is unknown, and PT-I is to be determined. Thus at time T-l, we need P~X~l, 1- -ID~-l'H~-l) "-~ a and ~e~-l is known.

279

Thus, the determination of optimum

max ~-l

PT-I

can be expressed as

~T-I

where YT-I : IT-I + E (yT IHT.I) with IT-! ^

In evaluating

E(YTI~T_ll , we use (7)' to express

8T

in terms of

ST_1

and ~T-I (PT-I) " Thus YT-1 is a function of ST-1 and PT-I ' and the indicated maximization can be carried out either analytically or approximately numerically. We sketch some approximation techniques for the ease where an analytical solution is not available in Section 5. General Case Defining the maxim~ value as price at time

k

YT~I ' etc., the determination of the optimum

can be expressed as

where

5.

DISCUSSIONS

The process of determining the optimal sequence of prices described in the previous section is generally rather complex, unless the marketeer chooses a simple criterion function such as setting to last stage, generally.

E(YTIHT_I)

Pt

by

E(xtlHt) = 0 . For example, in the next

cannot be expressed in an analytically closed form,

The situation is exactly analogous to that of Example C of Sec. III.3 of

Ref [3].

~ , which multiples

Pt ' appears as an unknown gain in the equation for

St

f

is

when

is linear and

J

given by

(2').

One way to simplify the price determination is to separate the price deterndnation from the estimation; for example, (certainty equivalent pricing policy), or p(elHk)

at time

k . Then

other words, at time

Pk

8

is replaced by 8

appearing in

is datez~Z_ned using

E(elH±) J

at time

is estimated using

8k as the parametem.

k , the demand function is approximately taken to be

xt(Pt) : f(Pt'St ) + ~t "

t In

280

Another approximation is based on the idea of predictive control [II]. for example that at time

t

the price

p

Assume

is maintained throughout the remainder of

the planning period (static price expectation).

This enables the specialist to have

his estimated sequence of stochastic excess demands for the remainder of the periods f(P'eT ) + ~T ' T = t,t+l, ..., T . Expand

k(.)

excess demand value up to quadratic terms, say.

E

+

and

,(.)

about this predicted

This gives an approximate value to

(+)

*(%)I s=t

as a function of

p . (This is a purely stochastic problem. )

Repeat the above for two other fixed

p

values, say

p + Ap

and

p - Ap

with some

Ap . Pass a q~dratic curve through these three points to obtain a quadratic

approximation to (+) as

~(p)

. Then perform

Max ~ ( p ) P The maximizing

p

is chosen to be

Pt " Repeat the whole process at

t+l .

In some dynamic systems, this approximation is known to give better results than open-loop feedback approximations.

See for example Sec. VII of [3] and [9,11] .

The details of vamious pricing policies such as certainty equivalent or openloop feedback pricing policies are left to a separate paper, where an intuitively appealing price adjustment mechanism can be shown to be derivable from (sub)optimization of certain criterion functions.

See [4].

We now briefly indicate modifications that are necessary when we drop the synchronized trading assumption. Let

R

be the basic period that the trading specialist uses in revising his

price, i.e., after each

R

period, he evaluates the past history of the excess

demands and decides to revise his price or not. Arrivals of buyers (and sellers) can be modeled by a stochastic point process, fore example by a Poisson process with intensity

p . The market demand then is

modeled by a compound Poisson process. Durdng the period Let

d(p) + {

R , suppose that the trading specialist observes

be the demand schedule where

~

with finite vamiance~ independent of the point process. The market excess demand~ now written as

z ~ is then

% N2 z : Nld(P) - N2s(p) + [i:l ~i - [j=l ~j where

~'s

and

~'s

are assumed to be independent.

N

buyers.

is a mean zero random variable

281

Since E(N I) = E(N 2) = pT wehave

E( )/pT = x(p) , where

x(p) : d(p) - s(p)

: (%2 ÷ o22)/pT

is the "relYresentative" customer's excess demand for the

eon~nodity. Therefore, if the specialist merely wants to devise a pricing policy which achieves zero excess demand in some prebabilistic sense, he could use Pt+l : Pt + Izt where

X

is related to

k , small positive adjustment coefficient by

For sufficiently small

k

and when

d(p)

and

s(p)

incompletely specified coefficients, it can be shown that such that

E(x(p*)) = 0

with finite

aspect will not be pursued f u ~ e r

I = k/QT .

are linear in E(Pt)

p

with

converges to

Var(x(p*)) under suitable assumptions.

p* ~his

here.

Another possibility is to treat the change in price as a statistical hypothesis testing~ zero hypothesis being that

x(p) = 0

and the alternate hypothesis being

x(p) # 0 . This is especially easy to apply when normally distributed.

~'s

and

q's

are assumed to be

See for example Ferguson, Lehman and Wilkes.

When the specialist wishes to use the criterion function (2), then must be added to

Ht . The reeursion formula for

p(81H t)

(NI,N 2 )

requires only a minimal

amount of modification, since the point process generating

N1

noise processes

Thus the formula for

~

and

q

are independent by assumption.

and

N2

and the 8t

quite analogous to (7) is obtained for this more general case. Under the constant future price expectation, one can express the probability distribution for the future excess demands which are used to obtain quadratic approximation to E(JIH t)

E(JIH t)

as outlined in Section 4.

generates an a p p r o x ~ t i o n

Min~zation

to the optimal

of this approxin~te

Pt "

The details must be left to a separate paper~ again because of the lengthy development. ~EFERENCES i.

Alchian, A.A., "Information Costs, Pricing and Resource Unea~loyment," in Phelps et el, Micro-economic Foundations of Employment and Inflation Theory, 27-52, W.W. No~ton and Company, Inc. (1970).

2.

Arrow, K., "Towards a Theory of Price-Adjustment," in The Allocation of Economic Resources, Abramo~itz, M. et al (ed.), Stanford University ~ s s

(1959).

282 3.

Aoki, M., 0 p t ~ a t i o n (1967).

of -Static

Systems, Chapter I I I , Academic Press

4.

Aoki, M., "On Some Price Adjustment Schemes of a Marketeer," presented at the 2nd Stochastic ContrDl Conference, NBER Conference on the Computer in ~onomic and Social Research, Chicago (June 1973).

5.

Barro, R.J., "A Theory of Monopolistic Price ~hjust~m.nt," Rev. Econ. Stud. 34 (i), 17-26 (January 1972).

6.

Clower, R.W., "Oligopoly Theory: A Dynamic Approach," Proc. 34th Ann. Conf. Western Economic Association (1969).

7.

Clower, R.W., "Some Theory of an Ignorant Monopolist," E e o n ~ c 716 (December 1959).

8.

Clower, R.W., "Competition, Monopoly, and the Theory of Price," Pakistan Economic Journal, 219-226 (September 1955).

9.

Dreyfus, S.E., "Some Types of Optimal Control of Stochastic Systems," SIAM J. Control 2_, 120-134 (January 1964).

Journal, 705-

10.

Hadar, J., and Hillinger, C., "Imperfect Competition with Unknown Demand," Rev; Econ. Stud. 366, 519-525 (1969).

11.

Tse, E., Bar-Shalom, Y., and Meier, L., "Wide-Sense Adaptive Dual Control for Nonlinear Stochastic Systems," to appear in IEEE ?PanS. Aut. Control.

PROBLEMS OF OPTIMAL ECONOMIC INVESTMENTS WITH FINITE LIF~IME CAPITAL

Bernardo Niccletti Istituto Elettroteonicc Universit~ di Napoli Napoli, Italy.

Luigi Marian~ Istitute di Elettrotecnica e Elettronica Universit~ di Padova Padova, Italy.

I. INTRODUCTION The theory of o p t i ~ m economic dynamic investments has been extensively studied. However,apart from simple models, the finite lifetime of capital has not been considered. As a matter of fact, in many allocation problems, either in the context of macro- or micro-economics, the effects of some decisions may have only finite duration in time. Aim of this paper is to focus this problem and to present a method for dealing with the problem in the discrete-time case. Mathematically, when one considers finite lifetime L, each decision affects not only the state of the systems at the following stage but also at a number of stages later. The usual prescription for handling problems of this type is to increase the original state space by L additional variables: as a consequence, even for problems with relatively short delays, the dimensionality of the problem is i~ creased beyond the reach of present-day computing capabilities and the new system is not controllable. D. G. Luenberger suggested an approach for exploiting the inherent structure of this problems in order to get a low-dimensional algorithm: the "cyclic dynamic programming". In the present paper, a similar approach is used

to refor~late the original

problem. However, rather than aiming directly to a computer algorithm, a "quasi-ps riodic discrete maximum principle" is derived, which allows one to write general canonical optimam conditions. From these, it is also possible to derive computer al gor it ~ . In ~ecti~n 2 and 3, as an example of applications, two models for systems w i ~ finite lifetime investments, in the context of macro- and micro-economy respectiv~ ly, are discussed. In Section 4, the general problem is stated. Using a mathematical programming approach, necessary conditions for the optimality are derived in

284

Section 5 for the general nonlinear case. For systems with linear dynamics, concave inequality constraints and convex performance index, the conditions derived are both necessary and sufficient and can be formulated as an extended discrete maximum principle. This is done in Section 6. Finally, in Section 7 the method is applied for deriving the optimality conditions for the models discussed in Section 2 and 3.

2. A NACROECONOMIC MODEL

To show the basic ideas of this paper, the simplest model for the accumulation of capital in one-sector closed econo,~" will be considered. As it has been suggested, (J. Tinbergen and H. C. Bos), when one considers economic growth, it is desiderable to distinguish between the stock of equipment or capital goods and the stock of capital. The difference lies in the fact that an industrial machine remains a constant volume of equipment until it is scrapped, whereas its contribution to the capital stock decreases because of its depreciation. As a consequence, the structure of the model is different from the classical one. Using the variable: t

=increasing ordinator, indicating the discrete-time instant considered;

k t =capital stock, at time t; Yt =national income (net product), at time t; it =gross investment, at time t; ct =consumption, at time t; b t =capital goods or volume of equipment, at time t; L r

=finite lifetime; I =--~-=constant proportionate rate of depreciation,

(where kt, Yt' it' ct' and b t are expressed in the same unit), the follo~Ang relationships hold (I)

Yt = ct + it

income identity

(2)

kt+ I = k t + it - r bt

gross investment identity

(3)

bt+ I = bt

(4)

Yt = f(bt)

production function

(5)

ct>/0 , i t ~ O

and bo, i I , .... , i_L given.

+

it

_

it_L

capital goods identity

285

Dynamic economic models are considered here as a tool in planning development. Particularly, the viewpoint considered is that of a government which is in a position to control the economy completely by choosing the gross investments and to plan perfectly so as to maximize a given social welfare function over the time T T-I (6)

~

wt(ct)

.

T is the planning horizon: in a logical sense, the relevant period of time is the entire future. However, development plans are often chosen to maximize a welfare function over a finite horizon T, subject to terminal stocks constraints which represent a weight attached to the welfare of the generation beyond the horizon

(7)

bT~%.

B~,,A MICROECONOMIC MODEL

In many dynamic allocation problems for

a

single firm, the effect of some de-

cisions may have only finite duration in time: this is the case, for example, in the planning of production plants with a finite period of economical operation, the replacement of Faculty members in educational institutions, the schedule for cutting and planting trees in a given area, etc. Problems of this type may be

described by

a discrete-time model of the following type. Using the variables: u t = capacity built, at time t; c(u) = cost of building capacity u; a

m interest rate;

x t = total capacity operating at time t; d

t

= demand at time t;

T

= planning horizon;

L

= plant lifetime;

(where ut, xt, and d t are expressed in the same unity), the following relationship hold (8)

xt+ I = x t + ut

t = O, ..., L-I

286

(9)

xt+ I = xt + u t - u t _ L

t = L, ..., T-I

(1o) xt>d t (11)

Xo =

t = I, ..., T

o.

It is desired to minimize the costs

over the planning horizon T

T-I

(12) 7

1 (l-~a)t

t=o

c(ut)

By choosing the building policy.

4. STATEMENT OF THE PROBLE M

The models discussed in Section 2 and 3 can be reformulated in a general way as follows. It is given a dynamic system governed by a set of nonlinear difference equations (13)

xk+ I = fk (Xk' Uk' k)

k = O, .... , L-I

(14)

Xk. I = fk (Xk' Uk' Uk-L' k)

k = L, .... , T-I

where k is an increasing ordinator (the stage number), x k and u k respectively

the

the state and control variables, and L (the finite lifetime) is less than T (the t~ me horizon). L and T are positive integers. (At a slight increase in notational ccm plexity, the whole discussion can be trivially extended to the case where both x k and u k are vectors). The problem is to find a sequence of ~Uk)

(k = O, .... ,T-I),

subject to belong to the constraint set represented by (15)

ek (xk, Uk' k) ~ O

k = O, ... ,

in such a way to minimize the performance index T-I

(16) J = >

~ (xk,

u k, k)

&..__

k=O starting from a given initial condition x o. Following Luenberger's approach, first one define the integer n such that (17)

(n-l) L ~ T ~ n L

and the vectors (the prime' denotes transpose)

287

(18)

x~- (x k, xk+L, ..., ~+(,_I)L)

k = O, I, ... , T-(n-I)L

(19)

x~ : ("k' xk+L' ""' xk+(n-2)L)

k ffiT-(n-I)L+I, ... , L

(20)

U~ = (~k' ~k.~' ""' ~k+(n-1)L )

k = O, I, ..., T-(n-I)L-I

(2~)

U~ = (~k' ~k.L' "'" ~÷(~-2)L )

k = T-(n-I)L,..., L-I

w i t h (18) and (20) of dimension n and (19) and (21) of dimension n-1. With these position, equations (13) and (14) can be expressed by the equivalent set

(22) xk÷I=F k(x k,U k)

k ffiO, I, ..., I.-I

while the constraints of type (15)

are

replaced by

(23) Ek(x k,Uk)>o

k = O, 1, ..., I.-I

and the performance index tc be minimized L-I

(24)

o:~

ak(~,u

k)

k=O where

(25) Gk(~,U k):

n-1

~÷iL(~÷iL'~',÷iL)

k

ffiO, I,... ,T-(n-I)L-I

i:0 n-2

(26)

ak(xk, ~k ) :

~_

~k÷iL (~k+iL' uk+iL)

k ffiT-(n-I)L,

... , L-I

i=O In addition, there is the consistency condition, which is reflected in the boundary conditions

(2?)

xjffi Xo X L]

which c ~ p l e s

'x°

giv~

the initial and the final states. In the c ~ e

T = nL the b ~ n d a r y

con-

dition is (28)

x ° ]=X ° i 5 Xo XL

given.

xT ]

The function Fk, Gk, and E k are supposed

to be

continuous and differentiable.

In the case of infinite horizon the preceding vectors have an infinite hum-

288

ber of components; what follows may be extended to this o~se by using functional an~ lysls (C~aon et al.).

~. MA~m~mTIOAL moa~m~nca APP~,O,AC~

This problem can be reformulated as a mathematical programming problem, by defining the vectors Z'

=

(X o, X I, .... ,XL, U o, U I, ... , UL_ I)

F(Z)' = (Fo-X I, FI-X 2, ..., F~_cX L)

E(z), =

(~o,

~I'

"'"

~'b.1 )"

The optimization problem becomes that of minimizing the scalar function G(Z) subject to F(Z)=O and E(Z)~O and to (27). By direct application of mathematical programming theory, it is possible to derive the following theorem. Theorem I In order that ~Uk} furnishes an optimal solution to problem (22)÷(24), it is necessary that there exist a real valued number ~, vectors )~k and/~k (of the same dimension of Xk+ I and ~

~Fk

~Gk

aFk

respectively), such that the following equations hold

Oa k

~Ek

~Ek

k=0,..,5-1

(30) ~-~k ÷ (~-~k)'xk"(~Tk)'~k = o

k=0,.. ,555551

(31) ~+i

k=0,..,L-1

=

Fk(~,uk)

(32) ~k(Xk,Uk)> 0

k=O,.. ,L-I

(~3) Z~k(Xk,Uk) = 0

k=0,.. ,555551

(34)

/%:~o

(35)

(o,~) Xo ~ xT.

(36) ~-1 = (o,i), "~5555-I (37) ~>~o

k--O,..,L-I

289

(+8) (4, q+,x~,.',x'L_+,~%,4;,-'e~_+) ¢ o. Constraint (38) in the statement of this theorem may b e omitted, and the variable~ be set equal to I, if the constraint qualification holds: for a general discussion of the conditions under which they hold, see Canon et al. In order to obtain sufficient conditions for the problem under consideration, further conditions over Fk,Ek, and Gk must be imposed. Using mathematical pro&,ramming theory t it is possible to show that the following theorem holds. Theorem 2 Given F k linear in ~

and Uk, given the initial condition Xo, given the per-

formance index to be minimized (24), with Gk convex in X k and Uk, and given the constraints (23) with E k concave in &

and Uk, in order that ~Uk~furnishes an optima]

solution, it is necessary and sufficient that there exist a real valued number~ and vectors ik and/~k such that equations (29) to (+8) hold. 6.A ~JASI-PERIODIC DISCRETE MAXIMUM PRINCIPLE The results of theorem 2 can be posed in the form of a discrete maximum principle. Define the Hamiltonian

(39) ~(xk,uk, Xk,~) =~Gk(Xk,U k) + X~Fk(X~,Uk)

.

k~O,..,L-1

Assume that the conditions of theorem 2 are satisfied and that (v°,u°,A@ ,%-o~ is ~k k k'~k ''I j the optimum solution. Proceeding in a way similar to Pearson and Sridhar, it is possible to show that (30) may be substituted by

(4o) Hk (X~o,UO,~O, k k ~o) = min ~k(X~'Uk'~'~ @) with

(+I) ~ = luk= mkCx~,uk)~o } . In terms of Hamiltonian, (29), (30), and (31) m y ~H k

~E k

(42) ~,k_l= ~---~&- (~--~&)',"k

be rewritten as follows

k--O,..,t.-1

+,& (4+)

3T k -

(44)

&+l

o o

"

"

If ~ is zi.ear in Uk, eq.. (43) must be s~betitutea by (40). Therefore t it is possible to conclude that conditions (42), (43) or (40),

290

(32) to (38), and (44) represent a discrete maximum principle• The particularity of this case lies in the consistency conditions (27) and (36); which account for the n~ me of "quasi-periodic discrete maximum principle" for its similarity to the prinoipie which holds for periodic systems (Bittanti, Fronza and Guardabassi).

7" EXAMPLES 7.1. An Economic Growth Problem. As a first example of application, the macroeconomic model outlined in Section 2 is considered. The production function f(bt) and the welfare function w(ct) are supposed to be concave monotone increasing. We also assume, for simplicity, that T=NL, so that the time horizon is integrally divided by the finite lifetime. Let

(50) C~ = (Ck,Ok+L,..,ck+(~_1)L)

~-

, ~

(bk,bk+L,..,bk+(~_1)i.)

~

(ik, ik+L,..,ik.(~_llL)

I,-I

;

~O,Z,..,U--1

N-I

(52) z~' = (ik_co,..,o) ,ik_L ~iven

.

k:O,I,..,L-1

The problem ca be re-stated as follows

(53) Bk+I = Bk*

A

I k - I~

k = O , . . ,I--1

k=O, . . , I , - 1

(54) Ik>O

(55) v(Bk)

zk>I0

k = O t . , tI,-1

(56) (0,..,0, I) BL >/ bT

(57)

b°]=

(o,i) ~o " (i,o) ~L

or

BL

(58) Wk(Ok) = Wk(F(B k) - Ik) with

(59) A -

I

0

0

°

•

-I

I

0

•

.

1

0

.

e

•

•

0-1 °

•

0 O 0 0 •

•

.

0

e

•

0

0-1

1

U s i n g theorem 2, w i t h I k , ~ k , and ~ k v e c t o r s of dimension N, ~ v e c t o r of dimension N-1

291 to take into account condition (57),.~scalar and ~ the scalar multiplier to take into account the condition (56), it is possible to derive the following necessary and su~ fioient conditions (canonical equations)

~Z(Bk)

(6o1

(61) Bk+ I = Bk + A ~ (62)

~L-I - (0,0,..,0,

~Wk(Ck)

~F(Bk) )÷ (

k=O,..,L-1 k=O,.. ,L-I

I~ I)'~

= - (I,0)'~

;A_I : - ( 0 , I ) ' #

~Wk(Ck) (63) 0 = - ~ ( O-~--k ) + A'A k - ~ k +fk

k=O,.. ,L-I

(64) F(Bk) - Ik>O

k=O,.. ,L-1

(65) ~(F(~ k) - zk) = 0

k=O,.. ,L-I

(66) A ~ o

k=O,..,L-1

(67)

~k>o

k=-O,..,L-I

(68) ~ Ik = 0

k=O,..,L-1

(69) ~k 3 °

k=O,.., 11-I

(70) ( o , o , . . , o , 1 ) (71)

(72) (73)

~L - b~ ~ o

~((o,o,..,o,~)

~L - ~T ) : o

~ o (o,z) B° = (z,o) BL .

For the presence of the boundary condition on the final value (70)t the condition (36)

becomes now

(52).

For actually computing the optimal solution, it is possible to employ an algorithm of the following type. Let us consider first of all the case that no one of the constraints (64) and (67) are binding. As a first trial (74)

Bo (I)

(I) ~o (with~(1) =oil I~1)do given with k k of dimension 2 for k:0,..,T-L-I and of dimension I for k--T-L,..,I,-I;~k of di mension 2 for k=-0,..,T-L and of dimension I for k~T-l~1,..~I.-S. A computational algorithm, similar to the one discussed for the economic growth model, may be used for actually finding the optimal solution.

8. CONCLUSIONS The optimum economic dynamic investments problem has been investigated in the case of finite lifetime investments. The problem is reformulated in such a way that, by using mathematical programming, it is possible to derive general canonical optimum necessary conditions. From these, a computer algorithm can he derived. For linear systems with convex performance index and concave constraints, the conditions found are both necessary and sufficient and may be posed in the form of a discrete maximum principle. As examples of application,

a

macroeoonomio growth model and a microeoonomic

plant investment model are considered.

Acknowledgment This work was supported by the Consiglio Nazionale Ricerche, CNR, Roma, Italy.

REFERENCES Bittanti~S., FronzatG.~ and @aardabassi,G., Discrete Periodic Optimization, Internal

294

Report 72-I~, Is%. Elettrotecnica ed Elettronioa, Pol. Milaao (1972). Canon,M.D., Cullum, C.D., and Polak,E., Theor[ of Optimal Cqntrol and Mathematical Programming, Mo Graw Hill, New York (1970). Luenberger,D.L., Cyclic Dynamic Programming: a Procedure for Problems with Fixed Delay, Operations Research, 19, 1101-1110 (1971). Pearson,J.B., and Sridhar,R. A Discrete Optimal Control Problem, IEEE Trans. on Auto. Con%____~.,AC-11, 171-174 (1966). Tinbergen, J., and Bos,H.C. Mathematical Models of Economic Growth, Mc Graw Hill, New York (1962).

SOME ECONOMIC MODELS OF MARKETS by M . ~ H ~ M o g r i d 6 e Centre for Environmental Studies, London*

I.

Introduction

2.

The Labour Market

3.

The Household Income Distribution

4.

The Car M~rket

5.

The Market for Space - the Land Market

6.

Conclusion References

i.

Introduction In this paper I am dra~#ing together a theme from a whole series of papers which I have written over the last few years. Several other authors have also contributed work along these lines, and I will be incorporating such work in this exposition. Indeed, the models which I am using have a very long sad distinguished history, both in mathematics, and more specifically in economic systems. The gamma distribution, first extensively tabulated by Pearson in 1928 has been used in its collapsed form of the negative exponential distribution in an economic system at least as far back as 1892 by Bleicher, while the logarithmic normal distribution has the distincZion of having a book written about its manifold uses in 1957 by Aitehison & Brown. What justifies a new look at these models, in my view, is the wealth of data now becoming available about the prices in markets with the advent of computers in governmental systems. The three major economic systems that I am going to discuss have all undergone, or are undergoing, a remarkable change in data availability and accuracy. This has made possible much more sensitive evaluations of the use of these models. In the first place we have the labour market. With the advent in 1964 of the graduated pension scheme in the United Kingdom, the Department of Health and Social Security now has a complete set of annual earnings records of everyone on one computer, and has been taking a I% completely random sample of these records annually since then. Since 1970, the Department of Employment has also been taking a 1% sample of the same records and obtaining details from employers of all matters concerning the earnings of each such employee in a given week. This again is now done annually. This has given us entirely new precision in our data on the labour market. In the second place, the car market, all licencing records both of cars and of individuals arc~ currently being placed on one computer by the Department of the Environment. The data which will be available here will hopefully include *Now with the Greater London Council, Department of Planning & Transportation.

296

not only the current ownership but also the previous car owned, but at present these details are only voluntary and we do not y e t have any indication of response rates. The third market that I wish to discuss is that of space, or more specifically the spatial location of households, worker residences and worker workplaces within urban areas. Here the key advance is that of geo-coding, or coordinate referencing, whereby each zone, ward or district in the urban area can be located with respect to all others. While such work has been done by hand before, the advent of the computer has made it possible to handle the enormous matrices that one gets when connecting worker residences with workplaces, eg. in London we have 1000 x lOO0 matrices. In London, coordinate referencing systems have been developed from the 1966 Census, both by the Department of the Environment and the Greater London Council for Census areas, and by Blumenfeld of University College London for the Traffic areas of the 1962 London Trafflc Survey as applied to 1966 Census data. Let us see then what statements we these models to markets. 2.

are

able to make about the applicability of

The Labour Market In a number of papers, I have shown that:

: 5.~

v 4 : ~.3 1

D

: 30 ml/min/mmHg

cmH20/l/sec,

o~ V 3 est le volume de rinqage et V 4 - V 3 le volume de l'~chantillon

et T A L e

temps

d' apn~e. La figure 3 concerne la courbe Fech = 9 (TA) pour les valeurs types. On remarque une l~g&re inflexion de la courbe vers 15 secondes. Elle est & rapporter changement

automatique de sensibilit~ introduit dans le syst@me pour la valeur la

plus faible des fractions ~chantillonn~es. corroborent

au

l'expression

analytique de la fraction ~chantillonn~e.

La pression barom~trique,

la r~sistance des voles a~riennes ainsi que le

d~bit expiratoire n'ont que peu d'influence La capacit~ de diffusion, les 3 coefficients

Les r~sultats obtenus par ~ s divers essais

£ondamentaux,

sum l'~volution du syst~me.

le volume r~siduel et le volume inspir~ sont

car ils influencent

tr~s £ortement la pente de la

371

d@croissance ainsi que la valeur inltiale de la fraction ~chantillonn~e. La fraction inspir~e et le d~bit inspiratoire n'ont d'influence que sur la valeur initiale du syst~me. La place et le vol~mle de l'~chantillon pr~lev~ ont une importance plus grande qu'on ne pourrait s'y attendre. En effet, si leur influence sur la valeur de la fraction alv~olaire en fonction du temps eat nulle (saul par augmentation ~ventuelle du temps d'apn~e) leur influence sur la valeur de la fraction @chantillonn@e est tr~s grande aussi bien en ce qui concerne la valeur intitale que la pente de la d@croissance de Fech = f (TA). Ceci a d@j~ ~t~ soulign~ par diff~rents auteurs et pourrait expliquer tout au moins en pattie la disparit@ des r~sultats ob~enus

au

moyen de dill@rents modules.

CALCUL AUTCMATIQUE DE LA CAPACITE DE DIFFUSION PULMONAIRE

L'expression analytique de notre module ne peut @tre invers~e pour en extraire la capacit~ de diffusion ; pour cette raison nous nous sommes attaches trouver un mode de calcul automatique de cette capacit~ A partir de la sortie du module et des param~tres r@ellement mesur@s. Nous proposons 5 modes de calcul tous bas@s sur le principe de l'acctunulateur analogique constitu~ de deux unit~s "track and store" mont~es en s~rie et reboucl@es sur l'entr~e (fig.4) cet acct~nulateur permet d'entreprendre la recherche pas ~ pas du coefficient de diffusion. La premiere idle est d'envoyer des increments proportionnels A la di£f6rence entre la fraction expir~e mesur@e et la fraction ~chantillonn@e calcul~e par le module. Le type de calcul (type I) repr~sent~ sur la Figure 5 ayant dur~ 43 p@riodes, noua n'avons pas recherch~ la proportion optimale qui aurait donn~ le temps de calcul le plus court. La seconde idle est d'envoyer des increments constants. Nous choisisons 3 types d'incr~ments : 0,1 UM - O,01 UM - 0,001 UM pour qu'on obtlenne des r~sultats avec tune erreur absolue de 10 -4 ou que l'accumulateur £ournisse

dix lois la valeur

de la capacit@ de diF£usion. Les quatre calculs pr~sent@s ci-apr~s sont des variations sur lea co~nutations entre chaq~e pas de calcul : - calcul de type @ La commutation au pas inf~rieur se fait d~s que Fech F E et la commutation au plus petit pas si 95 % < F E