Decomposition on the Basis of Structure of Discrete Systems. Application to Hierarchical Approach of a Water Quality Control Problem

Copyright
DECOMPOSITION ON THE BASIS OF STRUCTURE OF DISCRETE SYSTEMS. APPLICATION TO HIERARCHICAL APPROACH OF A WATER QUALITY CONTROL PROBLEM

J.

L. Calvet* and A. Mano**l

*Laboratoire d'Automatique et d'Analyse des Systemes du C.N.R.S., 7, Avenue du Colonel-Roehe, 31400 Toulouse, France **Instituto Superior Teenieo, Centra de Analise e Proeessamento de Sinais, Av. Roviseo Paris, Lisboa, Portugal Abstract : This paper deals with hierarchical algorithms for digital ~omputer control. Its purpose is to exploit the serial par~llel type of interconnection structure corresponding to a wide class of large scale systems encountered in many cases of practical interest (transportation systems, hydraulic systems, phys ico-clE mical processes, ... ) Several decomposition schemes are discussed for the optimal control of discrete dynamical systems and emphasis is put on convergence analysis and numerical application. They are also compared to the most popular class of algorithms proposed for this kind of problems by Tamura (gradient type algorithms based on duality property in convex problems). Numerical comparisons are made in the context of a constrained water quality problem in rivers. A detailed description of this problem is provided with some comments on the modelling of the process and its discretization. The control corresponds to the effluent discharge from treatment plants. The state constraints are spacially distribltP.d along the river stretch considered. This application shows the advantage of the proposed sequential schemes which lead to improvement greater than two order of magnitude in computation efficiency. Keywords: Discrete systems ; optimal control; partitioning coordinatiou ; distributed control; water quality control. INTRODUCTION This paper is de'l.Oted to the numerical solution of optimizat10u problem on large scale disc~e dynamical systems. Mainly concerned with the problem of dimensionality, it deals with: - iterative (indirect) methods for solving large scale equations - hierarchical approach of large scale systems and is directed toward distributed computing. One 0 f the purposes is to exploit some interconnection structure to provide usefull communication schemes _in order to coordinate the decisions ~ocally elaborated by the subsystems,in a decentralized fashion. In the reference (Calvet, 1981), we have discussed the application of such a methodology to the quadratic optimization on linear continous (and distributed-lag discrete)dynamical systems,interconnected in a serial-parallel pattern. In the present paper, we ~~ill limit ourselves to the discrete case .for which, we will mainly emphasize partitioning and convergence analysis as numerical application. I.This research has been partially supported by IN[CTunder contract nOI19.79.82 and by INIC. 139

decomposition-

Iterative methods of decomposition : The studied decomposition methods apply to the discnte case the double concept of approximation and iteration and are usually called: iterative methods, successive approximations method, indirect method, ... The approximation can be carried out by means of a Taylor series expansion (quasilinearisation f.i), a polynomial expansion (Pade rational approximation f .i) as by any other ways providing search directions. This leads to the construction of an approximating (k+1)t· sequence1x s. x(k+1)

T[x(k)]

where k~O is the iteration index and x(o) is the initial guess of x, which converges toward a fixed point of the mapping T, i.e. toward a solution XX of

Semi iterative methods : This concept, introduced by Varga (1962), defines non statiorary methods x(k+l) = Tk+l[x(k)] with a view to improving the conditions of applicability or convergence. Among the most popular, let us mention :

J. L. Calvet and A. Mano

140

- ~~_~l!~E~~!i~g_~iE~~!i~~~_~~!~~~~

Sequential decomposition with symmetric alternated directions x. (k+I/2)=T.[x.(k+I/2),x (k)] i=I~N 1 ~ J 1 j i

x(k+I/2)= TI [x(k)] x

(k+l)

T2 [x(k +1/2)J

x.(k+I)=T[

- ~~~_g~~~!~li~ed_~i£~~E~~~~_~~!~od~ (k+l) = (k) + F D r(k) x x K+I k+1 where:E. k + IC IR

~

Xj

i Xj

(k+I/2),

~

i=~1

(k+I) ]

(k+I/2l

X.

,

(k+I/2)

1

D + is a non singular matrix k 1 r

(k)

expresses the residual vector at SS.

step k.

1

Different decomposition algorithms can thereby be constructed after partitioning to iterRtively approximate a solution x~, taking into account such miscellaneous features as improvement of convergence, adaptation to the computing tool (parallel,vector~omputers), use of distributed computers, ... Let us mention the basic schemes illustrated by Fig.I-3 Parallel decomposition (Jacobi type algorithms). (k+ I) T. [x.(k)] jFi; i,j =I, ... ,N x. ~

~

J

X.

(k+l)

x

1

e(k+ I) ;0l> . 1

Fig.3 Distributed communication and local memorisation in an alternating directions process. OPTIMIZATION OF DISCRETE DYNAMICAL SYSTEMS Let us consider the optimization problem of type LQS N-I I T T Min L, (x Q x +u R u) n=O 2 n+1 n+1 n+1 n n n u

n=O, ... ,N-I

n

(I)

. UEIRP n

,

Xo given s .t. xn+1 = ~n xn +ru n n n=O, ... ,N-I Qn semi-definite positive n=I, ... ,N Rn definite positive n=o, ... ,N-I

Fig.1 Centralized communication parallel process.

~n

a

Sequential decomposition (Gauss-Seidel type algorithms) x.

(k+ I)

~

T [(k+l) i Xj

''1

(k)] . . J<1 ; 1> i i=I~N

Before to derive decomposition- coordination schemes exploiting the stair-case interconnection structure of the overall system,let us write as follows the stationarity conditions associated with (I) : =x - ip x -r u =0 n n-I n-I n-l n-I

n=l. .. ,N

=p -Q x - 4 \ =0 n n n n n+1

n=I, . .. ,N-I

(D

Jp

n

(l)

lx (0

I

n

x

=p -Q x =0 N NN N T

cD

"Tu

C

SS

l

Fig.2 Decentralized and scheduled communication in a sequential process.

(2)

n

=R u +r p =0 nnnn+1

n=O, ... ,N-I

For the sake of simplicity ani to better depict the partitioning,we will limit our formulation to iterative methods as previously defined. This means that we do not focus here on interesting semi-iterative methods sugges ted by cOn'-'exi ty ass ump tions and based on the l1inimum Principle as well on the D'. lalit:, theory. Thus, di rec t i te ra tio;:s for searching ~, In, lh can be derived from the overall expression of the stationarity conditions put into the form : Az=b ; A€.lR2Nmx2Nm ; b€.R2Nm (3) Let us recall that the concept of direct iterations is associated with the Prediction

Decomposition on the Easis of Structure of Discrete Systems coordination principle (Mesarovic and co\vorkers, 1970). Several block-decompositions can therefore be derived after a splitting A=AO-Al of the matrix A partitioned as follows

141

~ o

where the submatrices Aii are square and non singular : , .

,.

N'

A..E:Fn 1xn 1; L n'i = 2 Nm 11 i=l Denoting AD=Block-Diag[Aii], ~ and ~ respectively the strictly upper and lower part of A, different algorithms can be constructed in a classical way to iteratively approximate the unique solution of (2): AOz(k+l)= Alz(k)+b

(4)

The case N'=2 appears of a great interest since it leads to a matrix ~ in a decomposed form.

J

l2 Considering A=[All A A21 Al~

for that case,and the reducible form of All due to stair case structure, it is obvious that the first

- Parallel decomposition: AO =\ ; AI;: - (A +\) U - Sequential decomposition: Al'E-A

U

- SYT!1l'1etric alternating direction method: All a-AV AI2~-A

L

Not concentrating on the details of these methods (Varga 1962, Youne 1971, ... ) nor on the po ssible exploitation of the computer network, we will limit ourselves to the necessary and sufficient condition of convergence of these algorithms.

b100k in

tu"

[~1 ~~oan

be dimtly invw

ted by sequentially solving in forward Sense the state equations, and that the second block corresponding to AIIT can be directly inverted by a backward sequential resolution of the cos tate equations (i=N _ 1). These two procedures being operated in a parallel or a sequential fashion,in the overall decomposition-coordination algorithm (4). Theorem 2 - The parallel or sequential decomposition are either both convergent or both divergent. The necessary and sufficient condition to converge is:

Theorem I (Lions and co-workers,1974). The algorithm ~4) converges toward the unique solution z of (3) when k.... .., if and only if

T Al2 AII - A2 )<1 while a sufficient condition may be given by : IIAI211 I/A 21 H
I)P(T)<.1

verge, the sequential decomposition is asymptotically twice fast as the parallel scheme.

where

pis

l

the spec tral ratius of the i tera-

tion matrix T _ AD-IAl 2) if the , iteration matrix has some eigenvalues A= 1, then:

rg~ll-TJ =

r[A ll -

rg

[11 _~

2

Taking into account the specific structure of the considered matrix A, two kinds of partitioning are of interest, based on different orderings of the stationarity conditions (2). Partioning I.

zT= [xI" .xu PI" .PN ]

Decompositions based on the partitions hereafter yeld a separation between state and costate equations. -I T Pirst, denote Sn ~rn R n r n ; n=O, ... ,N-I

Proof: The proof for N.S.condition is straightforward from theorem I. The sufficient condition derives from the fact that F [·J~H -lIfor any norm an:! that "AlII! =.11. due to its reducible form where the diagonal blocks are identity matrices. The other considerations on convergence are usual for two-cyclic orderings,what is precisely the case when N'=2 (Varga 1962). Let us underline the interest of the sufficient condition which states that convergence may be always achieved for any coupling in between subsystems by conve niently choosing the ponderations matrices Rand Q in Al2 . 1 n n an d A21 respectlve y.

J. L. Calvet and A. Mano

142

Let us also observe that this decomposition follows a feasible path and may be implemented "on-line". The algori thm corresponding to sequential decomposition will be numbered (!)in the subsequent numerical application. At last, it seems interesting to compare for the considered partition, this iterative method wi th the well known Dual decentralized method by Tamura (1975). Let us recall that the l atter involves in a parallel fashion, a semi-iterative method to improve the costate variables using steepest descent direction: (k) Pn with d k+1 > 0 and Dk+1 definite positive Pn

(k+l)

=Pn

(k) 01+1 Dk+I'o + ( I

...

AlgOrithmm: Sequential decomposition (i=l .... N") Algorithm 4 : " " ( i = N ' .... I) Algorithm 5 : symmetric alternating direction method. A more detailed description of these algorithms is given in (Calvet, 1981) where the feasibilit~f the algorithm 3 is emphasized. Note also that algori thm 3' s known as "Cos ta te predic tion method" (Papageorgiou and Schmidt,1980 ;

'[Al

rn.

-1

Tl

11 R fn . nl1 2/ N' =N9~=Block-Diag n:o, ... ,N-1 -Q n +1 As opposite to previous decompositions, t IS partition involves the invertion of the ~ matrix. However, its interest lies in the bYock-tridiagonal form of A which is a case of optimal ordering and is particularly well - conditionned for overrelaxation, see (Varga, 1962) .

J

T{

Different decompositions experimented in numerical simulation have been numbered as follows:

Singh and Titli, 1978)

while the variables x and u are predicted by direct iter~tions ~especti­ vely from the remaining stationarity conditionsf. =0 ander. =0. Un xn Such a method with conjugate gradient will be referred as algorithm~ Partitioning 2. z = xIPI ... xNPN

cost criterium, such that the algorithm will converge if the autonomous subsystems are asymptotically stable.

foxo

o

Considering parallel decomposition, it must be also underlined that Al is only function of the ~ matrices, which means that relaxation, performeg throught these matrices, is therefore well-suited for the case of weak- c~pling between subsystems, e.g whenq>n= E.1? rv E. smaD.. Let us finally note that this parallel decomposition constitutes an adaptation to discrete systems of the algorithm by Takahara (1970). \VATER QUALITY CONTROL APPLICATION I. Statement

n~o

partitions visualized above appears very conv enient for decomposition in that case. 1/ N'= 2r-; --=+ AD'" Block-Diag[ll) =

-n

This decomposition into subsystems of minimal dimension with respect to this partitioning, has the advantage to avoid the invertion of the matrix ~. Analyzing convergence property for parallel decomposition, the necessar~ and sufficient condi tion states p[T }fll1-A]<1 while an usefull sufficieni condition may be given by

n'

P [T J 1~ max I~

j'fl

i ~N'

It~jl=

(I.1=P= I =* 't'

n

=4>n ;

1

1

1

assumed convex. Final effluents are discharged at ,foints liE..Ealong a river stretch :.t',,[O LJund er constant flow rates, with DO (dissolved oxygen) contents :C and BOD i contents:U .-u .. 1

The aim of the water quality control problem is to minimize the sun of all treatment costs constrained t o keeping DO concentration above some minimum l evel ~J lVithin the river s tre tch.t' .

n 1

when dealing with scalar equations: A;

remove u.of the BOD load at a cost J.(u.)

1

{I
The decomposition approach previously described is now applied to a constrained water quality control problem in rivers. A set of M waste sources is assumed to produce effluents lVith (BaD) a biochimical oxygen demand concentration Ui(i: I, ... ,M). Treatment plants

The model used to describe the spatial profiles of DO and BaD due to the steady effluent discharges is presented in the appendix. ~ 2

max n=I, ... ,N

of the problem and assumptions

Q =q ; r R n n nn

-I

T

rI'\. = -r

2

O'n n

The latter condition shows that it always exists a convenient choice of the weights coefficients qn and rn in the

).

2. Optimization on the spatial discretized model. A spatial discretized version of the problem is formulated as follolVs:

Decomposition on the Basis of Structure of Discrete Systems N-I

L: n=O

Min u :nE:.JO

algorithm 0,where double precision variables were required in order to achieve convergence near the solution, and reach a stop criterium smaller than 10- 3 .

J (x ,u )

n

n

143

n

n

Values for the other parameters were J

=lJ x T[qln 0 Jx n- 21 n+1 0 q2n n+1

rn = r =10.

subject to the model equations (see appendix):

Xl =4 ; x :OQ ui:-oo 1

x n+ I=J. 'n x n +rn un +z n ; x 0 given n=O, ... ,N-I and boundary cons train ts on bo th s ta te am con trol: I xn ~ ~I n=1 , ... ,N ui.:s

ui~ui

1

1i

n- I

- , ... , N- I

PN-QNxN+\ = 0

An (x.1 - xn I)

u.=sat ~

r_ L

= 0 .~ '

rni T r

Pni+

. n1

u.

~

where sat[.]

'>

n'"

i~

0') =0 i f (x _xl

.J

, n

J

)<.0

n=I, . .. ,N

[.]

-

u.

if [.].(. u

~

!----

•

A'CuJ-SO)

u.

~

i

max 10

l '

A

I

,..

n

"- 10

)

10

J

R

12

3

20

J

20

3

40

10

n

JI

3. Numerical r es ults and cor.nnents.

The optimization problem w?s solved on a computer HP 1000 of CAP$ (Lisboa), using the five referred approaches. l-lhen needed, costate variables were initialized to 0, while state variables were 3uessed as the ; values obtained from control u = U / l I 2 U.= U.:i>l. A stop criterium based on the ~ ~ (k) (k) (k- I ) error C ;; Ip - p was defined

I

as :" N

~ tn

n=1

--

"I

(k) ~ 10- 5 , except for Tamura

"2

.J

")

".5

10 . 1

10 . 4

121 x

If))

20

52 . 9

89.2

50 .

n,;

10

JI

I I

8

60

34.9

68.1

"".1

112 x 10J

o

110 JIS

100

19 . 4

48.8

59.5

79 x 10J

lOO

19.4

41 . 2

fio . 1

61 x 10]

5

--

110

8

0)

' .3

r70 . 2

0

CV

1

5.3

2 .9

f----

140.1

15.5

lA . '

9

813.2

12 . 1

94 .3

1396.6

122 .5

134.4

23

140. 1

151 . 9

20.1

C

13 . 6

C'

19.'

0

CV

6.2

I

t-------

3

o

eo ..

A'

Let us note at this point that the previously mentionned property of feasibility in s ome decompositions is no lon3er ensured.

n

~J~~ --

OpH ..,) C'ns t

20

----

~(k'»o

n

1, Il z 11] •

~ A

n

r-Optf...'. Control

T,.bh I - Spr.dfic= det .. , Intution" And (''''U

(k')+P..(k')(x -x 1)(k'+I)L

n

N

C

et ('1 2"-0, I)

The iterative GC ilCl1C S previously defined have been cumb ined 'vi th an i tera ti ve adaptation of the Kuhn"Tucker parameters (conjugate gratient a130rithm applied to a ARROl-l-HUR1HCZ-UZAWA type method) ;

. .- ".. "."-.---- - .n • t ••

,.

(".aliI!!

> u ~.

if u ~ [.] i

'V~;

2

~-.-

;i=1 , .. , ,M

[.] ~

n

i u i : ..-

Six cases A,A' ,B,C,C',D were" copsidered with a different number M and positions li of waste sourc~and different dimension N.Moreover, the treatment in the third plant was cons trained to be ~ 50 in case A', and q2 was let to 0.1 in case C', in order to reauce the state contribution in the criterium cost. Table 1 shows the specific datas for each case, as well as the solutions (see also fig.4-5) and the optimal cost. Table 2 shows the computation time (CPU) (s] for each case and each method. A better computational efficiency ~observed with the sequential decomposit~on ~ as well as with the symmetric alternating method ([) • Table 3 shows for cases C and D the results in terms of the number of iterations.

i=I, •.. ,H

+An =0

; VI"\.

Xo =10. ; Xo =0.

The modificatio ns induced by these constraints in the statiomrity conditi.ons (2) are given in the following where the ~n n=I, ... ,N are the Kuhn-Tucker parameters associated with state constraints. Pn-Q n x n _;t Tp n+1

; qln =q2n =1.

--

>2000

16 . 1

41.16

23.6

." Table 2: - CPU c.olaputlltfon t l .

~

C

D

(i)

120

114



1826

®

1286


liS

Allorith.

118]

5 .6 12 . '

482

112

Tahle J - MUlliheE' of iter.tion.

J. L. Calvet and A. Mano

144

strength of couplings between subsystems, to the choice of weighting coefficients in the cost function •... Hence the interest of numerical experiments of various decompositions on significant problems.

20

.

15

" "./"

o

'"

"',

""',

.....'

10

REFERENCES

o

'" 5

,:

j "

o.~----~~------.~------~ 0 .0 2.5 5.0 7.5

____~

10.0

Drouin, M. (1981) Elaboration de nouvelles structures de commande des processus complexes. These d'Etat, Univ. Paris XI

Dislance Fig.4 Optimal State response (case A)

Lions, J.L.et Marchouk,G.I.Editeurs(1974) Sur les methodes numeriques en sciences physiques et economiques. Dunod Ed.

20

"

10

o

.

, ,,"" ,,

15

: \:

:,"

\:

I

:.

.,

,\;\

I \

\

:

\:, \'

Calvet, J.L. (1981) Decomposition sur la base de la structure serie-parallele de systemes dynamiques. Congres AFCET Automatique, Nantes (France) 27-29 Octobre 1981

\

I

\

I

:'.: \

\

I

I' I

\'

\ . : \!

\,

\

I

\

\1

I

\1

~\

\

:\ :

: \:

\

Mesarovic, M.D, Macko, D. and Takahara Y. (19 7 0). Theory of hierarchical multilevel systems. Academic Press.

I', \

:

\

\

:, , \

\

'.:

:

\!

\

:'~~~ DO DOD 10

15

20

25

30

35

40

D ls la n c: e Fig.S Optimal state re~ponse (case D)

CONCLUSION This paper shows that application of splitting methods issued from Numerical Analysis, for the Hierarchical optimization of discrete systems, may be successful even when standard properties (srmmetric positive definiteness of A; A- \.ith non-negative elements such as M-matrices, Stieljes matrices, ... ), are not met- Exploitation of structural properties as physical interconnection between subsystems, exploitation of the different types of variables (state, costate and control variables), ... constitute a usefuUway to perform in a decentralized fashion the optimization of interconnected systems. It is shown in this paper, how optimization of discrete systems yields well-suited partitionings for decomposition and coordination. Let us also refer at this point, to the methodology proposed by Drouin (1981) in order to improve the local exploitation of the information on the process. Concerning the numerical applications, it must be stressed that if somme algorithms appeared very performant, however results of comparisons must be analyzed with some care due to the alteration of convergence prop.e rties by the heu ris tic approach used to solve constrained problems. These applications also show that comparisons between different partitions and relaxed schemes remain an open problem, with respect to the specific type of problems, to the

Papageorgiou, M. and Schmidt, G.(1980). On the hierarchical solution of non lin ear optimal control problems. J. of Large Scale systems Theory and Applications, Vol.l, n 0 4. pp265-271 Rinaldi,S., Sancini-Sessa ,R. Stehrest H. and Tamura, H. (1979). M.odell in g and control of river quality . Mac Graw Hill. Singh,H.G and Titli, A. (1978). Systems : decomposition, optimisation and control. Pergamon Press Tamura,H.(1975) Decentralized optimization for distributed-lag models of discrete systems. Automatica, vol.ll-pp593-602. Varga, R.S (1962) Matrix iterative analysis . Prentice Hall Young, D.M.(1971) Iterative solution of large linear systems. Academic Press.

APPENDIX - THE WATER QUALITY MODEL The steady state spatial profiles of th e water quality variables : DO (Dissolved Oxygen) and BOO (Biochemical oxygen demand) along a river stretch are usually described (Rinaldi and co-workers, 1979) by the linear differential model : dx(l) = F(l) x(l) + u(l) + v(l) dl x(O)=x

O

; l~~=

lo

(A . I)

L}

\fuere the comporents of the state vector (x(l) ",[xI (1) x (l)JT are respectively the 2 DO and BOO concentrations [mg][ l r I , v(1) is a vector of uncontrollable known inputs, and u (1) expresses a vector of distributed inputs associated with the effluent discharges.

145

Decomposition on the Basis of Structure of Discrete Systems

In accordance with the prob l em description in the previous section, u(l) takes the form:

(A.2)

where O( . ) is the unitary impulse function'~i and x (1.) are respectively the dilution rate and th e §tate vector at the upstream side of the discharge points 1., H is the number of treatment plants . 1 The feasible BOD treatmen t u. acts as a pointwise control, whil e the 1 BOD concentration U. o f the ef fluents is assumed to be fixed ~

-k

-k

2

v (1)

-

-k - k l 2

k4 V

'1

For numerical applications , th e following typical val ues we r e used : V = 2.lkm][dall (mean stream velocity) k = 0.3 [ day] -I (deoxygenation rate) l k2= 0.6 [day] -I (reaeration rate) k3= 0. 1 [day] k = 0.2[mgJ[l] 4 Cs

=

Ui

=

[da y]

rI

-I

1

=

(distributed BOD l oad)

(oxy ge n saturation concentration)

Cl J- I

4.[mg][1] -I

C.

0(.

-I (s edime ntation rate) -I

I 0. [ mg] [ 1 I00. [ mg]

i=I, ... ,M

0.3

1

Spatial discretization was considered with N+I discrete points nf~{o,I, ... ,N } such that each discharge point nifJ(Q: i=I, ... M ~~{nl' .~.,nH},~~ rresponds tu a discrete POlTIt,

JVoccK.

loe.

After integrations of eq. (A . I) over the interval ~n=ln+l-ln' w~ere A(l) and v (1) are assumed constant

on

each ~

n

,and

denoting x(1 ) = x , v(1 )=v , F(l )=F n n n n n n for 1 E.[l

n

F =

e

n

1

~

n+

=i>n

whe re

x +r u +z n n n n

;p n .=

IJ, the solution.

n

x + n

is written into the form

n=O, ...• N-l

F A

(I-c() e n

r n =_"e.Fn t:.n [°1

n

1 i=l, .... M

~n

J °

V

F(1)",

°

n+1

k2 Cs

l

V

V

x

Note thatfD ~

° only

d't

at the discharge points

nE.Jfo corresponding to the contribution of the effluents.