ELSEVIER
Contents lists available at SciVerse ScienceDirect
Advances i n Water Resources
j o u r n a l t i o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / a d v w a t r e sDe Saint-Venant equations-based model assessment in model predictive control
of open channel flow
M. Xu ^•*, R.R. NegenbornP.J. van Overloop ^, N.C. van de Giesen ^
^Section of Water Resources Management, Faculty of Civil Engineering and Geosciences, Delft University of Teclmology, Stevinweg 1, 2628CN, Delft, Tlie Netherlands
''Section of Marine and Transport Teclmology, Faculty of Mechanical, Maritime and Materials Engineering, Delft University ofTeclviology, Mekelweg2,2628CD, Delft, The Netherlands
A R T I C L E I N F O A B S T R A C T
Article history: Received 11 July 2011
Received in revised form 3 July 2012 Accepted 5 July 2012
Available online 25 July 2012
Keywords:
Model predictive control Quadratic Programming Open channel flow
Sequential Quadratic Programming Saint-Venant equations
Model predictive control (MFC) is a model-based control technique that uses an optimization algorithm to generate optimal control acdons. Based on the model used in optimizadon, MFC approaches can be categorized as linear or nonlinear. Both classes have advantages and disadvantages in terms of control accuracy and computational time. A typical linear model in open channel water management is the Inte-grator Delay (ID) model, while a nonlinear model usually refers to the Saint-Venant equations. In eariier work, we proposed the use of linearized Saint-Venant equations for MFC, where the model is formulated in a linear time-varying format and time-varying parameters are estimated outside of the optimization. Quadratic Programming (OF) is used to solve the optimization problem. However, the control accuracy of such an MFC scheme is not clear. In this paper, we compare this approach w i t h an MFC scheme that uses Sequential (Quadratic Programming (SOP) to solve the optimization problem. Because the estimation of the time-varying parameters is integrated in the optimization in SQ.P, the solutions from SQF-based MFC are expected to be superior to the solutions of QF-based approach, However, SQP can be computa-tionally expensive. A simulation experiment illustrates that the QF-based MFC approach using a linear-ized Saint-Venant model has an accurate approximation of the control performance of SQP.
Crown Copyright © 2012 Published by Elsevier Ltd. All rights reserved.
1. Introduction
Over the last decade, m o d e l predictive c o n t r o l (MFC) o f open channel f l o w has been a subject o f extensive study [ 1 - 6 ] . MFC is a model-based c o n t r o l technique t h a t uses an o p t i m i z a t i o n algo-r i t h m to genealgo-rate o p t i m a l c o n t algo-r o l actions. Advantages o f MFC aalgo-re t h a t i t predicts the f u t u r e system dynamics, therefore being able to take i n t o account f u t u r e k n o w n disturbances. I t can also deal w i t h constraints w i t h i n the o p t i m i z a t i o n . Based on the type o f m o d e l used i n the o p t i m i z a t i o n , MFC approaches can be catego-rized as linear or nonlinear. The focus o f MPC i n open channel w a t e r management is m a i n l y o n e f f i c i e n t w a t e r delivery i n i r r i g a -t i o n sys-tems, and river opera-tions f o r flood or d r o u g h -t p r e v e n -t i o n . A c o m m o n feature o f the exisring research is that, typically, linear models are used f o r p r e d i c t i n g the system dynamics, such as the reservoir m o d e l and the classical Integrator Delay (ID) m o d e l i n [ 1 - 5 ] . Under certain assumptions, these linear models can approx-i m a t e the nonlapprox-inear system dynamapprox-ics w e l l . The MPC o p t approx-i m approx-i z a t approx-i o n problems w h e n using such linear models are easy and fast to solve. Moreover, guaranteed global o p t i m a l solutions can be f o u n d .
* Corresponding author. Tel.: +31 0 152782345; fax: +31 0 152785559. E-mail address: min.xu@tudelft.nl (M. Xu).
A n o n l i n e a r m o d e l can n o r m a l l y include m o r e system dynamics than a linear one. This extra i n f o r m a t i o n i n the n o n l i n e a r m o d e l may increase the c o n t r o l accuracy i n MPC. However, due to the use o f such a n o n l i n e a r model, the o p t i m i z a t i o n p r o b l e m can be-come non-convex and hard t o solve. Indeed, this is the case w h e n using the Saint-Venant equations. Theoretically, a guarantee f o r f i n d i n g the global o p t i m u m f o r nonlinear o p t i m i z a t i o n cannot o f -ten be g i v e n [ 7 ] . Since the o p t i m a l action needs t o be taken w i t h i n a prescribed t i m e p e r i o d i n r e a l - t i m e c o n t r o l , c o m p u t a t i o n a l t i m e is i m p o r t a n t i n achieving the o p t i m u m . U n f o r t u n a t e l y , such a n o n -linear MPC scheme can be very time consuming, e.g., due t o t h e CPU-intensive m o d e l executions f o r the n u m e r i c a l calculation o f gradients o f a Lagrangian f u n c t i o n w i t h respect t o the c o n t r o l v a r i -ables, especially i n t h e areas w h e r e these gradients are flat. This c o m p u t a t i o n a l c o m p l e x i t y i n MPC using such a n o n l i n e a r m o d e l was also stated b y Barjas Blanco [ 6 ] . Therefore, t h e y used a series of reservoir models instead.
Some researchers use a d j o i n t sensitivity analysis t o speed u p the n o n l i n e a r o p t i m i z a t i o n b y analytically c a l c u l a t i n g t h e g r a d i -ents o f the Lagrangian f u n c t i o n w i t h respect to the c o n t r o l variables [8,9]. This is attractive f o r m a k i n g such an MPC i m p l e -m e n t a t i o n feasible i n r e a l - t i -m e c o n t r o l , b u t i t needs extensive analytical analysis o f the nonlinear m o d e l a n d its derivatives beforehand. Moreover, any change to the c o n t r o l p r o b l e m requires a n e w analytical d e r i v a t i o n . For these reasons, the
0309-1708/$ - see front matter Crown Copyright © 2012 Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10,1016/j.advwatres.2012.07.004
38 M. Xu et al./Advances in Water Resources 49 (2012) 37-45
a d j o i n t s e n s i t i v i t y analysis is n o t conducted i n this paper. I n -stead, i n [ 1 0 ] , w e proposed an MPC scheme u s i n g linearized Saint-Venant equations i n a t i m e - v a r y i n g f o r m a t to a p p r o x i m a t e t h e n o n l i n e a r dynamics. The scheme requires a m u c h m o r e c o m -plex d i s c r e t i z a t i o n and m a t h e m a t i c a l f o r m u l a t i o n o f t h e equa-tions t h a n t h e reservoir m o d e l i n [ 6 ] . This MPC scheme solves t h e o p t i m i z a t i o n p r o b l e m w i t h a standard Quadratic Program-m i n g (Q.P) solver, w h i c h considers the Program-m o d e l constraints as linear.
The MPC approach i n [10] is f o u n d t o be the m o s t accurate f o r c o m p a r i s o n w i t h MPC u s i n g an Integrator Delay m o d e l and a Re-duced Saint-Venant m o d e l . The proposed m e t h o d f o r m u l a t e s the Saint-Venant equations as a linear t i m e - v a r y i n g state-space m o d e l . I t uses a ' F o r w a r d E s t i m a t i o n ' t o estimate the t i m e v a r y i n g p a r a m -eters outside o f the o p t i m i z a t i o n , based o n t h e o p t i m a l solutions over a p r e d i c t i o n h o r i z o n f r o m the previous c o n t r o l step. However, due t o t h e lack o f i n f o r m a t i o n at the last p r e d i c t i o n step, the o p t i -m a l solutions i n the previous c o n t r o l step are n o t o p t i -m a l any-more i n t h e present step. Therefore, i t is unclear w h a t the performance o f this QP-based MPC c o n t r o l l e r is. The purpose o f t h i s w o r k is t o ex-plore the accuracy o f t h e c o n t r o l procedure o f [10] b y c o m p a r i n g the results w i t h an MPC c o n t r o l l e r t h a t f o r m u l a t e s Sequential Quad r a t i c P r o g r a m m i n g (SQP) problems anQuad solves the entire t i m e -v a r y i n g Saint-Venant equations w i t h i n the o p t i m i z a t i o n . A c c o r d i n g to S c h i t t k o w s k i [ 1 1 ] , SQP is a state-of-the-art m e t h o d f o r s o l v i n g n o n l i n e a r p r o g r a m m i n g problems. Here the MPC scheme using this m e t h o d is called SQP-based MPC.
I n t h i s paper, w e focus o n the p e r f o r m a n c e assessment o f the t w o MPC schemes i n t e r m s o f w a t e r level deviations f r o m t h e t a r -get and t h e c o n t r o l actions. Because o f t h e i n t e g r a t e d calculation o f t h e t i m e - v a r y i n g parameters w i t h i n the o p t i m i z a t i o n , the solutions f r o m SQP-based MPC are expected to be superior to the solutions of QP-based approach, g i v e n s u f f i c i e n t c o m p u t a t i o n t i m e . I t is the question h o w the t w o m e t h o d s compare i n t e r m s o f c o m p u t a t i o n a l t i m e and c o n t r o l accuracy. Regarding the c o n t r o l accuracy, the SQP-based MPC can be used as a benchmark. I n a d d i t i o n , f o r the QP-based MPC, iterations are added b e t w e e n the ' F o r w a r d Estimat i o n ' and Estimathe 'QuadraEstimatic P r o g r a m m i n g ' blocks, i n order Estimat o c o m p e n -sate t h e i n f l u e n c e o f the lack o f i n f o r m a t i o n at the last p r e d i c t i o n step. Therefore, another goal o f this w o r k is t o investigate the sig-nificance o f t h i s i n f l u e n c e .
This paper is organized as f o l l o w s . Section 2 describes the m a i n components o f MPC, i n c l u d i n g the open channel flow m o d e l i n g and t h e o p t i m i z a t i o n p r o b l e m f o r m u l a t i o n . I t summarizes the QP-based MPC scheme u s i n g the linearized Saint-Venant m o d e l and introduces t h e SQP-based MPC scheme. Section 3 introduces the test case used to compare the c o n t r o l p e r f o r m a n c e b e t w e e n the t w o MPC schemes. A detailed d e m o n s t r a t i o n o f the results is g i v e n i n Section 4 and conclusions and f u t u r e research are g i v e n i n Section 5.
2. Model predictive control of open channel flow
M o d e l p r e d i c t i v e c o n t r o l has a general structure w h i c h uses an i n t e r n a l m o d e l t o p r e d i c t the f u t u r e system dynamics over a finite p r e d i c t i o n h o r i z o n and solves a constrained o p t i m i z a t i o n p r o b l e m w i t h a certain o p t i m i z a t i o n a l g o r i t h m . MPC uses online o p t i m i z a -t i o n , w h i c h means -the o p -t i m i z a -t i o n is conduc-ted a-t every c o n -t r o l t i m e step and o n l y the first c o n t r o l action over the p r e d i c t i o n h o r i -z o n is a p p l i e d t o the system. A t y p i c a l MPC c o n t r o l p r o b l e m i n open channel w a t e r m a n a g e m e n t is t o m a i n t a i n a w a t e r level dis-t a n dis-t d o w n s dis-t r e a m o f a c o n dis-t r o l sdis-trucdis-ture adis-t dis-t h e end o f dis-the canal reach, w h i c h w i l l also be the subject o f this paper. I n the f o l l o w i n g sections, w e discuss the m a i n components o f MPC f o r such a sys-t e m : i n sys-t e r n a l m o d e l and o p sys-t i m i z a sys-t i o n .
2.3. Open channel flow model
I n order t o c o n t r o l the open channel flow w i t h MPC, the d y n a m -ics o f the system needs t o be p r o p e r l y d e f i n e d i n the i n t e r n a l m o d e l of the controller. Open channel flow dynamics is usually described by the Saint-Venant equations, w h i c h c o n t a i n the mass and m o m e n t u m conservations [12] s h o w n i n Eqs. (1) and ( 2 ) :
dA„
9 Qdd d((lv)
dt dx
'dx
•g m = 00)
(2)w h e r e A„ is the w e t t e d area [ m ^ ] , Q is the flow [m^/s], qi is the lat-eral i n f l o w per u n i t l e n g t h [ m ^ / s / m ] , v is the average flow velocity [ m / s ] , w h i c h equals QJA^. g is the w a t e r level above the reference plane [ m ] , Q is the Chezy coefficient [m^'-^/s], R is the hydraulic ra-dius [ m ] , w h i c h equals AJPf{Pf is the w e t t e d perimeter [ m ] ) and g is the gravity acceleration [m/s^], A t is t i m e step [s] and A x is spa-tial increment [ m ] .
According to Stelling and D u i n m e i j e r [ 1 3 ] , the Saint-Venant equations can be spatially discretized w i t h staggered grids. A s e m i - i m p l i c i t scheme is applied to the t i m e integration, w h e r e the advection t e r m i n the m o m e n t u m equation is e x p l i c i t l y discretized by a firstorder u p w i n d m e t h o d . The f r i c t i o n t e r m is l i n e -arized b y using |Q| e x p l i c i t l y . A l l other t e r m s are i m p l i c i t . I n this way, the Saint-Venant equations are linearized at every t i m e step. S u b s t i t u t i n g the velocities o f step n + 1 f r o m the discretized version of Eq. (2) i n t o Eq. (1), the w a t e r levels can be calculated w i t h a t r i -diagonal system, and the velocities are updated w i t h the calculated w a t e r levels t h r o u g h the m o m e n t u m Eq. (2). The detailed discret-i z a t discret-i o n o f the Sadiscret-int-Venant equatdiscret-ions discret-is discret-i n c l u d e d the A p p e n d discret-i x .
In general, there is no specific f o r m a t f o r the m o d e l constraints i n MPC. However, i n QP-based MPC, the i n t e r n a l m o d e l is usually f o r m u l a t e d as a linear state-space system. Due t o the inter-connec-t i o n b e inter-connec-t w e e n w a inter-connec-t e r levels and velociinter-connec-ties, inter-connec-the Saininter-connec-t-Venaninter-connec-t equa-tions are a p p r o x i m a t e d b y a linear state-space m o d e l t h a t is t i m e - v a r y i n g as s h o w n i n Eq. ( 3 ) :
yk+l k„k Qkjk
( 3 ) w h e r e x is the state vector, A is the state m a t r i x , u is the control i n -p u t vector, Bu is the c o n t r o l i n -p u t m a t r i x , d is the disturbance vector, Bd is the disturbance m a t r i x and k is the t i m e step index. The de-tailed f o r m u l a t i o n o f each m a t r i x is i n c l u d e d i n the Appendix.
2.2. Genetic MPC formulation
Typically, an MPC p r o b l e m i n open channel w a t e r management solves the m i n i m i z a t i o n o f a quadratic objective f u n c t i o n , subject t o linear or n o n l i n e a r m o d e l e q u a l i t y constraints and linear i n e q u a l i t y constraints o n t h e c o n t r o l i n p u t s . The reason to use a quadratic objective f u n c t i o n is t o balance b o t h positive and nega-t i v e varianega-tions o f snega-tanega-tes and c o n nega-t r o l inpunega-ts, such as w a nega-t e r level deviations f r o m the target level and the change o f c o n t r o l l e d struc-t u r e flow. The f o r m u l a struc-t i o n can struc-t h e n be w r i struc-t struc-t e n f o r a cerstruc-tain c o n struc-t r o l t i m e step k: (-n-l I J - O s.t.h,(X\ [ƒ") = 0 1 = 1, riiX^U") s^O ! = 1 , . , xt+j+i
•^[("'""•'ywuu'+j]
j-o J (4) w h e r e J represents a quadratic objective f u n c t i o n , X' ..., x"*"]^ and U" = [u" u"*"-^]^ are the states and control inputs over the p r e d i c t i o n h o r i z o n w i t h a length o f n, h, and r,- are the i t hequality and i n e q u a l i t y constraints, respectively, trie is the number of equality constraints, m, is the n u m b e r o f i n e q u a l i t y constraints,
Wx and Wu are diagonal matrices representing the w e i g h t i n g factors
on the state x and the control i n p u t u, respectively.
I n open channel f l o w system, t y p i c a l l y the state x''*-**' = ^ is the w a t e r level d e v i a t i o n f r o m the target level ic,t) and the c o n t r o l i n p u t u^+J = - QJ+^"' is the change o f structure f l o w , w h e r e Qc is the c o n t r o l l e d s t r u c t u r e f l o w . The d o w n s t r e a m w a t e r level o f a canal reach is assumed c o n t r o l l e d , and Wx o n l y penalizes the c o n t r o l l e d w a t e r level d e v i a t i o n . The equality con-straints reflect the dynamics o f the system, e.g. the Saint-Venant equations i n this case. The structure f l o w Qj- is restricted t y p i c a l l y between the m i n i m u m and m a x i m u m f l o w s o f Qc.min and Qr.max as the i n e q u a l i t y constraints, n a m e l y Qc,min < Qf^"* =^ Qc.max, 0 ' = 0 , . . . , n - l ) . Note t h a t b o t h e q u a l i t y and i n e q u a l i t y constraints are over the p r e d i c t i o n h o r i z o n n.
A c c o r d i n g to the t y p e o f m o d e l constraints, the obtained o p t i m i z a t i o n p r o b l e m can be solved w i t h d i f f e r e n t solvers. For e x a m -ple, i f the m o d e l constraints are linear, Quadratic P r o g r a m m i n g (QP) can be use. I f the m o d e l constraints are nonlinear. Sequential Quadratic P r o g r a m m i n g (SQP) can be used. I n this paper, Eq. ( 3 ) is a linear t i m e - v a r y i n g system w i t h proper discretization. It is a linear a p p r o x i m a t i o n o f partial d i f f e r e n t i a l equations at each t i m e step w i t h t i m e - v a r y i n g parameters. However, the e q u a t i o n is actually nonlinear over a f i n i t e horizon, w h i c h exactly needs t o be consid-ered i n MPC as the i n t e r n a l m o d e l . I n the f o l l o w i n g sections, w e discuss t w o m e t h o d s to estimate the t i m e - v a r y i n g parameters i n MPC.
2.3. QP-based model predictive control
QP-based MPC is i n t e n d e d to solve a quadratic o p t i m i z a t i o n p r o b l e m subject to linear constraints. X u et al. [ 1 0 ] approximates the linear t i m e - v a r y i n g Saint-Venant m o d e l over the p r e d i c t i o n h o r i z o n as a linear m o d e l . Fig. 1 shows the w o r k f l o w o f the QP-based MPC c o n t r o l l i n g a w a t e r system. I n this MPC scheme, t h e t i m e - v a r y i n g parameters o f A'', and are calculated outside of the o p t i m i z a t i o n solver and estimated t h r o u g h a ' F o r w a r d Esti-m a t i o n ' procedure. ' F o r w a r d E s t i Esti-m a t i o n ' siEsti-mulates the Saint-Ve-n a Saint-Ve-n t equatioSaint-Ve-ns over a f i Saint-Ve-n i t e p r e d i c t i o Saint-Ve-n h o r i z o Saint-Ve-n based o Saint-Ve-n the previous o p t i m a l solutions w h i c h , i n this paper, are the o p t i m a l gate flows over the p r e d i c t i o n h o r i z o n .
However, due to the c o n t r o l step change f r o m k - 1 to k, the o p t i m a l flow o f the last p r e d i c t i o n step is missing. Therefore, the o p t i m a l solutions at step - 1 are n o t o p t i m a l any m o r e at step k w i t h the n e w system i n f o r m a t i o n o f t h e last p r e d i c t i o n step. This lack o f system i n f o r m a t i o n is severe f o r the v e r y first c o n t r o l step w h e n no previous o p t i m a l solutions are available at a l l . I n order to handle this d r a w b a c k and increase the c o n t r o l accuracy, i t e r a -tions b e t w e e n the blocks o f ' F o r w a r d E s t i m a t i o n ' and 'Quadratic P r o g r a m m i n g ' can be i n c l u d e d i n the procedure. I n practice, n o t all the w a t e r levels are measured along a canal. Therefore, a stan-dard K a l m a n filter [ 1 4 ] is used to estimate the u n m e a s u r e d w a t e r levels w i t h a l i m i t e d n u m b e r o f measurements.
W h e n the m o d e l Eq. ( 3 ) over a p r e d i c t i o n h o r i z o n is s u b s t i t u t e d i n t o t h e objective f u n c t i o n ( 4 ) w i t h calculated A'', Bl and B^, the o p t i m i z a t i o n p r o b l e m can be w r i t t e n as:
Full water levels and velocities as initial condition
K a l m a n I'ilier W a i c r S y s t e m
Measured water levels
\ /
Optimized control Iblw of the first prediction step
i-onvard E s t i m a t i o n L i n e a r inequality constraints No
L
Full water levels and velocities as inilial condition
S a i n t - V c i i a n t M o d e l L J
Optimized control How over the prediction
horizon
O b j e c t i v e (unction
O P S o l v e r
Full water levels and velocities over the prediction horizon L i n e a r limc-varvin!.! S a i i i t - V e n a i U m o d e l Yes Q u a d r a t i c I Prognimiiniit; I O P - h a s e d M P C
40 M. Xu et ai/Advances in Water Resources 49 (2012) 37-45
min/(U") =
minr^(t/YH"i;''-O
(5)
n - l
w h e r e H is the Hessian m a t r i x , ƒ is the Jacobian m a t r i x and gc is a constant m a t r i x .
Because the constrained m o d e l is linear w i t h p r e - d e t e r m i n e d and f i x e d t i m e v a r y i n g parameters o f A'', and B^, the Hessian and the Jacobian can be a n a l y t i c a l l y calculated w i t h one m o d e l r u n over t h e p r e d i c t i o n h o r i z o n [ 3 ] . Therefore, there is no need f o r extra m o d e l executions f o r the g r a d i e n t c o m p a r i n g t o n u m e r i -cal -calculation, w h i c h saves a large a m o u n t o f c o m p u t a t i o n a l t i m e , especially w h e n i t e r a t i o n s are used w i t h i n the o p t i m i z a t i o n . But the MPC procedure does need (niter + 1 ) x 2n m o d e l executions f o r the ' F o r w a r d E s t i m a t i o n ' t o estimate the t i m e v a r y i n g parame-ters and f o r the c a l c u l a t i o n o f the Hessian and Jacobian, w h e r e Hjter is t h e n u m b e r o f i t e r a t i o n s b e t w e e n the ' F o r w a r d E s t i m a t i o n ' and ' ( i u a d r a t i c Programming". The o p t i m i z a t i o n p r o b l e m o f Eq. (5) is a t y p i c a l (Quadratic P r o g r a m m i n g p r o b l e m , w h i c h is i n our case solved b y the standard MATLAB f u n c t i o n 'quadprog' [ 1 5 ] .
2.4. SQP-based model predictive control
I n SQPbased MPC approach, a specific m o d e l f o r m a t is n o t r e -q u i r e d . However, because o f the linear t i m e - v a r y i n g state-space m o d e l setup f o r the QP-based MPC, t h i s m o d e l setup is also a p p l i e d i n a SQP-based MPC scheme. Fig. 2 illustrates the w o r k f l o w o f an SQP-based MPC c o n t r o l l i n g the same w a t e r system as i n Section 2.3. Compared t o QP-based MPC, the 'QP-based M P C block i n Fig. 1 is replaced b y 'SQP-based M P C . I n this scheme, the entire m o d e l over a f i n i t e p r e d i c t i o n h o r i z o n is integrated i n a sequence o f QP sub-problems, w h i c h are solved b y a QP o p t i m i z a t i o n solver u n t i l the objective f u n c t i o n value converges. Because the t i m e -v a r y i n g parameters are calculated i n t e r n a l l y i n the o p t i m i z a t i o n , t h e solutions are guaranteed t o be o p t i m a l .
I n Sequential Quadratic P r o g r a m m i n g (SQP), the objective f u n c t i o n is r e q u i r e d t o have a second order d e r i v a t i v e and the
constraints are required t o be c o n t i n u o u s l y d i f f e r e n t i a b l e . SQP f o r -mulates a sequence of QP sub-problems based o n the Hessian a p p r o x i m a t i o n o f the Lagrangian f u n c t i o n b y t a k i n g the second order Taylor expansion of the Lagrangian f u n c t i o n , and uses an i t e r -ative w a y t o solve the QP sub-problems. I n this paper, the MATLAB f u n c t i o n ' f m i n c o n ' [ 1 5 1 8 ] is used to solver t h e SQP p r o b l e m . H o w -ever, f o r the sake o f f a c i l i t a t i o n o f e x p l a i n i n g the n u m b e r o f m o d e l executions r e q u i r e d ( w h i c h is related t o the c o m p u t a t i o n a l t i m e ) , w e p r o v i d e some relevant equations o f SQP.
SQP is based o n the Lagrangian f u n c t i o n s h o w n i n Eq. (6) w i t h the Lagrange m u l t i p l i e r s o n the m o d e l constraints.
L(X' , [ / \ A ' , N ' ' ) = ; ( X \ t / " ) -w h e r e L is Ilk _ r , , k + l ( 6 ) and the Lagrangian f u n c t i o n , = [X( , . . . , 1 Nf = , . . . , ii^^"]^ are the Lagrange m u l t i p l i e r s .
The QP sub-problems at c o n t r o l step k are f o r m u l a t e d according to Eq. (7) b y t a k i n g the second order Taylor expansion o f the Lagrangian f u n c t i o n . The solutions o f the QP sub-problems S are used to f o r m u l a t e the next i t e r a t i o n u n t i l the objective f u n c t i o n values converge.
min UX X « , U ' , A * , N «
mm s"
s.t. fl,(Xp, Up) -H Vxfii(Xp, Lip)s, + Vu/i,(Xp, lJp)S2 = 0, i = 1, ri(Xp,Up)-FVxr,(Xp,Up)s, -H V(;r,(Xp,[/p)S2 sc 0, 1 = 1 , . .
(7) , JTli w h e r e H^ ^ is the Hessian m a t r i x o f t h e Lagrangian f u n c t i o n at
) . ƒ / ' . is the Jacobian m a t r i x o f the Lagrangian f u n c t i o n at " p ) . &.P is t h e constant m a t r i x , S''^ is
t a i n i n g s\ = X f
-X;, and Si
" — ' ^" - u L , - ui Note t h a t the subscript p represents the i t e r a t i o n index i n SQP.Because o f the d i f f i c u l t y i n a n a l y t i c a l l y calculating the gradient of the Lagrangian f u n c t i o n , a n u m e r i c a l m e t h o d is o f t e n used t o estimate the g r a d i e n t by executing the p r e d i c t i o n m o d e l . F o r w a r d d i f f e r e n c e is one o f the n u m e r i c a l m e t h o d s f o r the e s t i m a t i o n by assuming a s m a l l p e r t u r b a t i o n 5 on the state or the c o n t r o l i n p u t , as described i n Eq. ( 8 ) :
Wnier Sv^ieni
Mea.surcd waier levels
Optimized control tlow of tiic first prediclion step
Ralnian Kilter
Full water levels and velocities as initial condition L i n e a r liiiic-varyitig Saini-Veiiniii model O b i c c l i v e fiint'iion W Q P suh-prohlein L i n e a r itiet|iiality const mints
^
Q P S o b e r No Inilial guess S y P b a s e d M F C-Fig. 2. SQP-based MFC controlling a water system.
v * I l k A' M*!
V * l i t M ' i
r (8)
The Lagrange m u l t i p l i e r s i n calculating the gradient o f the Lagrangian f u n c t i o n can be obtained t h r o u g h the K a r u s h K u h n -Tucker (KICF) conditions s h o w n i n Eq. (9). KICT conditions are t h e first order necessary o p t i m a l i t y conditions t h a t need t o be f u l f i l l e d f o r nonlinear o p t i m i z a t i o n [19,20]. I t is obtained b y t a k i n g the par-t i a l d e r i v a par-t i v e o f par-the Lagrangian f u n c par-t i o n w i par-t h respecpar-t par-t o X'', U'', A!' and N".
VL(X^,U^,Aj,Nj) = Vx/(X^\[/^)
4 k• V t J ( X ^ , U j )
A lik\ h ^ A t V x h , - ( X ^ , Ul) + J2K^ijr>iXl, Ul) = 0
i=l i=l
A!;,h,(x^,u^) = o
(9) N: pi - p / " > 0I n order to n u m e r i c a l l y calculate t h e Jacobian m a t r i x o f t h e Lagrangian f u n c t i o n , n(2n + 1 ) m o d e l r u n s are required. These m o d e l runs i n c l u d e n runs t o calculate t h e Lagrangian f u n c t i o n v a l -ues o f L^ML^k+n and L„t+n-i and 2 x n x n a d d i t i o n a l runs to calculate the Lagrangian f u n c t i o n values w i t h small p e r t u r b a t i o n on each o f the variables o f state X and c o n t r o l i n p u t U, n a m e l y L^k+T+s,... ,Lxktn+sLuk^s,... ,Luk+n-i+i. Because a l l the m o d e l variables over the p r e d i c t i o n h o r i z o n i n X and U are dependent, a s m a l l
p e r t u r b a t i o n on one variable w i l l change the m o d e l constraints and influence all Lagrangian f u n c t i o n values. Thus, the model needs to be executed over the p r e d i c t i o n h o r i z o n n f o r each small p e r t u r b a -t i o n o f -the 2n variables i n order -to calcula-te -the Jacobian m a -t r i x . For the Hessian m a t r i x , Quasi-Newton m e t h o d is used to calculate the Hessian a p p r o x i m a t i o n . The B r o y d e n F l e t c h e r G o l d f a r b -Shanno (BFGS) m e t h o d is one o f the m o s t p o p u l a r algorithms to update the Hessian a p p r o x i m a t i o n based on the Jacobian m a t r i x [ 1 9 ] . Therefore, there is no need f o r extra model executions to cal-culate the Hessian m a t r i x o f the Lagrangian f u n c t i o n .
In total, n(2n + 1) m o d e l runs are required f o r each o f the itera-tions i n the n o n l i n e a r o p t i m i z a t i o n w i t h the f o r w a r d difference m e t h o d .
3. Test case
In order to assess the performance o f the QP-based MPC scheme using the linearized Saint-Venant m o d e l i n comparison w i t h the SQ.P-based MPC scheme, w e w i l l p e r f o r m s i m u l a t i o n experiments. Fig. 3 describes a test canal reach and includes all the geometric parameter values. I n the test case, w e assume no lateral f l o w s and t r y to c o n t r o l the d o w n s t r e a m w a t e r level o f the canal reach at - 3 . 2 m MSL ( m e a n sea level) by operating the d o w n s t r e a m gate w i t h a m a x i m u m f l o w of 4 m^/s, given an upstream i n f l o w d i s t u r -bance. Fig. 4 shows the i n f l o w disturbance scenario, i n v o l v i n g large and f r e q u e n t variations. The discharges are above the m a x i m u m d o w n s t r e a m gate f l o w i n certain periods.
T w o s i m u l a t i o n experiments are setup to compare the QP-based and SQP-based MPC schemes. Figs. 5 and 6 s h o w the procedures. I n e x p e r i m e n t (a) (Fig. 5), b o t h MPC schemes use the same i n f o r m a -t i o n f r o m -the w a -t e r sys-tem, and o n l y -the o p -t i m a l c o n -t r o l ac-tions f r o m SQP-based MPC are sent back to the w a t e r system. This e x p e r i m e n t illustrates the controllers' behavior w i t h d i f f e r e n t pred i c t i o n mopredels i n terms o f the calculation o f timevarying p a r a m -eters. I n e x p e r i m e n t ( b ) (Fig. 6), b o t h QP-based MPC and SQP-based MPC f o r m t h e i r o w n loops w i t h the w a t e r systems. The t w o w a t e r systems start f r o m the same i n i t i a l c o n d i t i o n . This e x p e r i m e n t is i n t e n d e d to show the influence o f the controllers on the t o t a l w a t e r system loop. I n a d d i t i o n , according to the QP-based MPC scheme s h o w n i n Fig. 1, d i f f e r e n t numbers o f iterations are tested i n each MPC i m p l e m e n t a t i o n , especially f o r the first c o n t r o l step w h e n no previous c o n t r o l i n f o r m a t i o n is available f o r calculating the t i m e - v a r y i n g parameters o f the state-space model.
Each e x p e r i m e n t covers a s i m u l a t i o n period o f 20 h w i t h a s i m -u l a t i o n t i m e step o f 2 s. The gate is operated every 4 m i n -using the o p t i m a l c o n t r o l action. The reach is discretized i n 500 segments, r e s u l t i n g i n 500 w a t e r level states. A Kalman filter is i m p l e m e n t e d to estimate the unmeasured w a t e r levels and velocities as the i n i -t i a l c o n d i -t i o n o f -t h e m o d e l cons-train-ts based o n -the measured w a t e r levels at the upstream and d o w n s t r e a m side o f the reach (details o f the Kalman filter are i n [10]). The o p t i m i z a t i o n p r o b l e m in this case is o n l y subject to the Saint-Venant equations and there
400 eOO 800 1000 1200 Time (minute)
Fig. 4. Upstream inflow disturbance scenario.
E.spcrimenl (u)
\V:ilei .Svsleni SOI'-based MI't
Ql"-bas(.'d MPC
Fig. 5. Experiment of controller behavior.
Experiment (b)
i
Water System SQP-based MPC
Water System SQP-based MPC
Water System QP-based MPC
Water System QP-based MPC
Fig. 6. Experiment of control influence on the water system.
is only one structure i n c o n t r o l . Therefore, the n u m b e r o f e q u a l i t y constraints (me) and the n u m b e r o f i n e q u a l i t y constraints (m,) b o t h equal t o 1.
4. Results
In t h i s section, w e p e r f o r m the experiments described i n t h e previous section to assess the performance o f t h e QP-based and SQP-based MPC schemes. I n a d d i t i o n , w e w i l l p r o v i d e an analysis of the c o m p u t a t i o n a l t i m e requirements.
.4m" means the elevation is 1.4m txjlmv mean sea level
42 M. Xu et al/Advances in Water Resources 49 (2012) 37-45
4.1. Results of control performance
4.1.1. Experiment of controller behavior
Before s t a r t i n g e x p e r i m e n t (a), i t is necessary t o analyze the convergence o f the o b j e c t i v e f u n c t i o n values o f the QP-based MPC scheme u s i n g the linearized SaintVenant m o d e l , since i t e r a -tions b e t w e e n t h e ' F o r w a r d E s t i m a t i o n ' and 'QP solver' are used to increase the m o d e l accuracy. Over the iterations, the p e r f o r m a n c e is expected t o increase, as the u n k n o w n i n f o r m a t i o n f r o m the pre-vious o p t i m i z a t i o n is updated. Fig. 7 shows t h i s convergence o f the objective f u n c t i o n values at c o n t r o l steps 1, 13, 67, 115 and 193, w h e r e u p s t r e a m i n f l o w s are above the m a x i m u m d o w n s t r e a m c o n t r o l f l o w a n d the c o n t r o l l i m i t s m a y be reached. A t these c o n -t r o l s-teps, -t h e u p s -t r e a m i n f l o w s i n -the p r e d i c -t i o n also change sig-n i f i c a sig-n t l y , as is i l l u s t r a t e d i sig-n Fig. 4. Note t h a t the f i r s t objective f u n c t i o n values i n Fig. 7 are w i t h o u t iterations. I t is clear t h a t the objective f u n c t i o n values converge a f t e r a certain n u m b e r o f itera-tions. B u t iterations o n l y have significant i n f l u e n c e f o r t h e very first c o n t r o l step, w h i l e the influence on the rest o f the c o n t r o l steps is negligible. Fig. 8 also shows the e v o l u t i o n o f predicted
w a t e r levels and gate flows o f QP-based MPC over 30 iterations at the first c o n t r o l step. The e v o l u t i o n results are s i m i l a r t o Fig. 7 and s h o w a fast convergence.
A direct i n d i c a t o r to compare the behavior o f the controllers is t o analyze the c o n t r o l goal o f m i n i m i z i n g the objective f u n c t i o n values. Fig. 9 shows t h e results o f e x p e r i m e n t (a). Because o f t h e lack o f previous o p t i m a l i n f o r m a t i o n i n the QP-based MPC scheme, especially at the v e r y first c o n t r o l step, the figure also plots the re-sult o f this MPC scheme w i t h 10 iterations between the ' F o r w a r d Estimation' and 'QP solver' at the first c o n t r o l step. A m a x i m u m of 10 i n i t i a l iterations is taken as tests showed t h a t t a k i n g a larger n u m b e r o f iterations does n o t alter the solutions any f u r t h e r (Fig. 7). Fig. 9 shows t h a t SQP-based MPC has r e l a t i v e l y l o w costs per c o n t r o l step over the p r e d i c t i o n h o r i z o n and o u t p e r f o r m s the QP-based MPC i n t h a t sense. The difference between the objective f u n c t i o n values o f t h e t w o MPC schemes shows t h a t the lack o f i n f o r m a t i o n f r o m the previous c o n t r o l step does influence the c o n -t r o l performance, especially f o r -the very firs-t c o n -t r o l s-tep, oscilla-tions occur due t o the complete ignorance o f t h e previous i n f o r m a t i o n at the b e g i n n i n g . B u t the influence vanishes q u i c k l y
Number of iterations
Fig. 7. Convergence of objective function values of the QP-based MPC scheme at
control steps 1, 13, 67,115, and 193.
SQP-based MPC QP-based MPC: no iteration QP-based MPC: 10 iterations 50 100 150 20D
Control step (4 minutes)
Fig. 9. Comparison of objective function values between Sbased MPC and
QP-based MPC (10 iterations are only applied at the first control step).
15 30
^
25 in nee 10 '8 iteratio n differ e%
15 •B ag e*
l o l •cen i 5 5 CD Q. + 5 10 15 20 25 Prediction horizon (step / 4 minutes)100 150 200 Control step (4 minutes)
Fig. 8. Evolution of the predicted water levels and discharges of QP-based MPC over
30 iterations at the first control step.
Fig. 10. Percentage difference in objective function values between SQP-based MPC
as the s i m u l a t i o n proceeds. On the other hand, the QP-based MPC can f o l l o w the SQP-based MPC performance t r a j e c t o r y very w e l l .
Moreover, there is a m a x i m u m o f o n l y 14.39% difference i n the objective f u n c t i o n values b e t w e e n SQP-based MPC and QP-based MPC w i t h 10 i n i t i a l iterations, s h o w n i n Fig. 10. The percentage is calculated by t a k i n g the difference o f their objective f u n c t i o n values d i v i d e d by t h e objective f u n c t i o n value o f the SQP-based MPC scheme. The larger differences m a i n l y happen at the begin-n i begin-n g o f the s i m u l a t i o begin-n , w h e r e m a begin-n y f l o w cobegin-nstraibegin-nts are activated, b u t no previous i n f o r m a t i o n is available as w i l l be s h o w n i n Fig. 12 later on. A l t h o u g h t h e c o n t r o l f l o w also reaches the constraints i n the s i m u l a t i o n afterwards, f o r example at about 800 m i n s h o w n i n Fig. 12, the percentage difference between the t w o MPC schemes is small (<5%) due to the previously available i n f o r m a t i o n (at about c o n t r o l step 200).
4.1.2. Experiment of control influence on the water system A f t e r assessing the behavior o f the controllers i n e x p e r i m e n t (a), i t is necessary to analyze t h e i r c o n t r o l influence on the w a t e r sys-tems over a f u l l s i m u l a t i o n . Regarding the control performance, the
1 1.5 2 2.5 Simulation step (2 seconds)
3.5 x i o "
Fig. 11. Controlled water levels in SQ.P-based MPC and QP-based MPC without
iterations. 5 4.5 4 3,5 g 2.5 O 2 1,5 - r 1-• S Q P - b a s e d M P C - Q P - b a s e d M P C U p s t r e a m flow > li ii 1 Ii 1 1 - , i. I (i i i '
n
l l l p '\,'r I \ M'I \i m
I iV
' I .' / \ r 10
0,5 1 1,5 2 2,5 3 3,5Simulation Step (2 seconds) ^ I D ' '
Fig. 12. Controlled discharge in SQP-based MPC and QP-based MPC without
iterations (upstream inflow is plotted to indicate the downstream flow trends).
f i r s t indicators are the c o n t r o l l e d w a t e r level and discharge e v o l u -t i o n over -the s i m u l a -t i o n horizon, s h o w n i n Figs. 11 and 12 f o r b o -t h MPC schemes, respectively. Both controllers can c o n t r o l the w a t e r levels very close to the target, and the system behavior using these t w o controllers is very similar. The w a t e r level oscillations t h a t are visible i n Fig. 11 are realistic and caused by the discrete steps t h a t the controller takes every 4 m i n , w h i l e the f l o w dynamics is s i m u -lated using a 2 s t i m e i n t e r v a l . The m a x i m u m c o n t r o l f l o w o f 4 m^/s is reached several rimes d u r i n g the w h o l e s i m u l a t i o n as s h o w n i n Fig. 12.
The performance o f the c o n t r o l l e d w a t e r levels and discharges can also be indicated t h r o u g h q u a n t i t a t i v e performance indicators, such as the M a x i m u m Absolute Value ( M A E ) i n d i c a t i n g the per-centage o f m a x i m u m w a t e r level d e v i a t i o n f r o m the target and the Integrated Absolute Discharge Change (lAQ) calculated by the integral o f the absolute f l o w changes over the s i m u l a t i o n m i n u s the absolute f l o w difference o f the first and last s i m u l a t i o n steps [21 ] . l A Q reflects the tear and wear o f the gate over the w h o l e s i m ulation. The comparison o f the t w o MPC schemes using these i n d i -cators is s h o w n i n Table 1. As can be seen, t h e differences i n c o n t r o l performance b e t w e e n t h e schemes are s m a l l . Due t o the higher m o d e l accuracy, SQP-based MPC has s l i g h t i y l o w e r values f o r the M A E and l A Q indicators, w h i c h indicates s l i g h t i y better performance.
4.2. Results of computational time
C o m p u t a t i o n a l t i m e is i m p o r t a n t f o r r e a l t i m e c o n t r o l i m p l e m e n t a t i o n . Before the n e x t c o n t r o l time step is reached, the o p t i -m i z a t i o n results should be available. Table 2 shows the c o m p u t a t i o n a l time o f d i f f e r e n t c o m p o n e n t s i n b o t h MPC execu-tions. According to the MPC procedures i n Figs. 1 and 2, the ' t o t a l t i m e per c o n t r o l step' represents the average c o m p u t a t i o n a l t i m e f r o m ' w a t e r system' to ' w a t e r system' t h r o u g h ' K a l m a n filter' and 'MPC block. The 'control t i m e ' is the specific time f o r the c o n t r o l generation i n ' M P C block. The ' m o d e l calls i n c o n t r o l ' represents the t i m e o f m o d e l simulations w i t h i n the ' M P C block. For QP-based MPC, the ' m o d e l calls i n c o n t r o l ' includes t h e time o f c a l l i n g the m o d e l i n the 'Forward E s t i m a t i o n ' and the controller, w h i l e f o r SQPbased MPC i t is the time o f c a l l i n g the m o d e l i n t h e SQP o p t i -m i z a t i o n . The calculation of b o t h MPC sche-mes was o n a 6 4 - b i t c o m p u t e r w i t h an 8 Gb i n t e r n a l m e m o r y .
From Table 2, i t f o l l o w s t h a t 99.37% o f t h e t o t a l c o m p u t a t i o n a l time i n SQP-based MPC is used to generate t h e c o n t r o l actions.
Table 1
Performance indicators of both SQP-based MPC and QP-based MPC. Maximum absolute value
(MAE") (%)
Integrated absolute discharge change (lAQ) (m^/s)
SQP-based -0.55 13.52
MPC
QP-based ^0.60 13.55
MPC
MAE is negative because the target level is negative.
Table 2
Computational time components in both QP-based MPC and SQP-based MPC executions.
SQP-based MPC QP-based MPC Total time per control step (s) 1761.25 54.39
Control time (s) 1750.14 44.40
44 M. Xu et al/Advances in Water Resources 49 (2012) 37-45
Table 3
Number of model executions per control step in QP-based MPC and SQP-based MPC.
QP-based MPC SQP-based MPC
No Five iteration iterations
n,reriterations
Number of 2n (5 + l ) x 2 n (n,fer + 1) X 2n n{2n -i-1) x üueT model
executions
w h i l e 99.9% o f t h e c o n t r o l actions calculation t i m e is spent on the m o d e l executions w i t h i n the o p t i m i z a t i o n . This illustrates the i m p o r t a n c e f o r r e d u c i n g t h e t o t a l n u m b e r o f m o d e l executions and speeding u p the m o d e l c o m p u t a t i o n s . I t is noted t h a t the ' c o n -t r o l -t i m e ' and -t h e ' m o d e l calls i n c o n -t r o l ' do n o -t include -the -t i m e o f K a l m a n f i l t e r i n g .
Moreover, there exist large c o m p u t a t i o n a l t i m e differences be-t w e e n be-t h e be-t w o MPC execube-tions. Since be-t h e acbe-tual c o n be-t r o l be-t i m e sbe-tep is 4 m i n u t e s , t h e table indicates t h a t the o p t i m i z a t i o n t i m e o f SQP-based MPC is unacceptable i n this case, w h i l e the i m p l e m e n t a t i o n of QPbased MPC is possible and 2 or 3 iterations can even be a l -l o w e d . The t i m e o f ' m o d e -l ca-l-ls i n c o n t r o -l ' is n o t the m a j o r t i m e c o n s u m p t i o n i n QP-based MPC and m o s t o f the c o m p u t a t i o n a l t i m e is spent o n the large m a t r i x m u l t i p l i c a t i o n s i n order t o b u i l d up the Hessian and Jacobian over the p r e d i c t i o n , according t o X u e t a l . [ 1 0 ] .
Because m o d e l c a l l i n g takes m o s t o f the c o m p u t a t i o n a l t i m e i n SQP-based MPC, i t is i n t e r e s t i n g t o analyze the n u m b e r o f m o d e l executions. Table 3 compares t h e n u m b e r o f m o d e l executions per c o n t r o l step o f the t w o MPC schemes, w h i c h is a f u n c t i o n o f the n u m b e r o f i t e r a t i o n s Hjter and the p r e d i c t i o n h o r i z o n n. For the QP-based MPC, niter is t h e n u m b e r o f iterations b e t w e e n the 'Forward E s t i m a t i o n ' a n d 'Quadratic P r o g r a m m i n g ' blocks. For SQP-based MPC, t h e n u m b e r o f iterations per c o n t r o l step is taken as an average n u m b e r o f iterations over the f u l l c o n t r o l steps (fiiter). It should be noticed t h a t t h e m o d e l execution i n QP-based MPC o n l y occurs i n t h e ' F o r w a r d E s t i m a t i o n ' and preparation o f Hessian and Jacobian matrices b u i l d u p over the p r e d i c t i o n h o r i z o n . There-fore, the n u m b e r o f executions increases l i n e a r l y w i t h t h e n u m b e r o f iterations b e t w e e n t h e 'Forward E s t i m a t i o n ' and 'Quadratic Pro-g r a m m i n Pro-g ' . Because t h e Hessian and Jacobian matrices are pre-cal-culated, there is no need t o execute the m o d e l w i t h i n the QP o p t i m i z a t i o n , a l t h o u g h QP i t e r a t i v e l y solves the o p t i m i z a t i o n p r o b
-70
10 I ' ' ' ' '
0 50 100 150 200 250 300
Control step (4 minutes)
Fig. 13. Number of iterations used in SQP optimization.
l e m . Fig. 13 shows t h e n u m b e r o f iterations used i n t h e SQP i n SQP-based MPC. There are fiiter = 45.86 iterations per c o n t r o l step on average, and every i t e r a t i o n requires n ( 2 n + 1) m o d e l calls t o calcu-late the Hessian a p p r o x i m a t i o n and the Jacobian matrices. For t h a t reason, the c o m p u t a t i o n a l difference b e t w e e n the t w o MPC schemes is significant.
5. Conclusions and future research
I n this paper, a comparison has been made between QP-based and SQP-based MPC schemes f o r c o n t r o l o f open channel flow t h a t can be described by 1-dimensional Saint-Venant equations. T w o experiments w e r e conducted f o r testing the behaviors o f c o n t r o l -lers and t h e i r influence on the w a t e r systems. Both MPC schemes presented i n t h i s paper can c o n t r o l the w a t e r system very w e l l . Due t o the integrated calculation o f the t i m e - v a r y i n g parameters o f t h e Saint-Venant equations, more accurate p r e d i c t i o n m o d e l is achieved i n SQP-based MPC. Furthermore, SQP-based MPC achieves better c o n t r o l performance regarding the m i n i m i z a t i o n o f the objective f u n c t i o n . On t h e other hand, SQP-based MPC is m o r e c o m p u t a t i o n a l l y expensive. The QPbased MPC using the l i n -earized Saint-Venant equations was tested i n a single canal reach i n this paper. However, i n general, as the c o m p l e x i t y o f the p r o b -l e m increases, the c o m p u t a t i o n a -l t i m e w i -l -l increase according-ly w h e n the same m o d e l accuracy (spatial discretization) is reached. In QP-based MPC using linearized Saint-Venant model, itera-tions b e t w e e n the 'Forward E s t i m a t i o n ' and 'Quadratic Program-m i n g ' can i Program-m p r o v e the c o n t r o l results. But the i Program-m p r o v e Program-m e n t is o n l y s i g n i f i c a n t at the very first c o n t r o l step, w h e n there is no pre-vious o p t i m a l i n f o r m a t i o n at all and m a n y constraints are active. The i n f l u e n c e o f t h i s lack o f previous o p t i m a l i n f o r m a t i o n dies out q u i c k l y as the s i m u l a t i o n proceeds and the c o n t r o l results converge fast over t h e iterations. As a b e n c h m a r k f o r c o n t r o l p e r f o r -mance, the results o f SQP-based MPC s h o w t h a t the procedure o f QP-based MPC w i t h the 'Forward E s t i m a t i o n ' is an effective and e f f i c i e n t w a y t o deal w i t h n o n l i n e a r i t y o f the m o d e l constraints.
This paper also provides an i n t e r e s t i n g finding t h a t m o d e l exe-cutions take 99.9% o f the c o n t r o l t i m e i n the o p t i m i z a t i o n o f SQP-based MPC. This suggests a w o r k i n g d i r e c t i o n f o r SQP-SQP-based MPC i n the f u t u r e regarding r e d u c i n g the t o t a l n u m b e r o f m o d e l execu-tions and speeding u p the m o d e l calculaexecu-tions, f o r example, the ad-j o i n t m e t h o d . I n a d d i t i o n , s w i t c h i n g all the calculations to l o w level p r o g r a m m i n g w i l l also decrease the c o m p u t a t i o n a l t i m e s i g n i f i -cantly, since t h e m a t r i x inverse i n the m o d e l calculation takes m u c h t i m e i n MATLAB calculation. Because o f the n o n - c o n v e x i t y of the o p t i m i z a t i o n , there is no guarantee o f global o p t i m u m i n nonlinear o p t i m i z a t i o n . Once the c o m p u t a t i o n a l b u r d e n o f SQP-based MPC o p t i m i z a t i o n is conquered, a m u l t i - s t a r t m e t h o d , a n u m b e r o f SQP-based MPC executions w i t h m u l t i p l e i n i t i a l values, can be i m p l e m e n t e d t o increase the chance o f reaching the global o p t i m u m .
M o r e practically, another f u t u r e research d i r e c t i o n regarding MPC is to analyze t h e i n f l u e n c e o f f u t u r e u n c e r t a i n t y on the p r o -posed QP-based MPC procedure. For example, w h a t is predicted i n t h e previous step does n o t happen or w h a t is n o t predicted i n t h e previous step does happen, etc. U n l i k e the theoretical w o r k i n this paper, w h i c h assumes perfect predictions, real w o r l d i m p l e -m e n t a t i o n s have t o deal w i t h these changing predictions.
Acknowledgement
This research was supported p a r t l y b y an I B M 2010-2011 Ph.D. Fellowship A w a r d and p a r t l y by the VENI project " I n t e l l i g e n t m u l tiagent c o n t r o l f o r flexible c o o r d i n a t i o n o f transport hubs" ( p r o -j e c t 11210) o f the D u t c h Technology Foundation STW.
Appendix A. Discretization of the Saint-Venant equations and
formulation as linear time-varying state-space model
Tlie deptli-averaged Saint-Venant equations are as f o l l o w s :
öK
dQ. dt ^ ö x •Qiö t + dx ^ ^ ' " d x ^ ^ d - R - A .
(A.2)
The equations can be discretized w i t h staggered grids i n space and s e m i - i m p l i c i t f o r t i m e integration as (the advection t e r m i n (A.4) is associated w i t h the positive f l o w ) :
-A, At "i+t/2 "1+1/2 Q-l-1/2 Q i + 1 / 2 Ax 1 ( A 3 ) A t w,i+1/2 Ax ''i+l/2 Ax
+g
J ( + l Ax - + g '^/+l/2l''/+l/2 l ^ f + l / 2 l ' " i + 1 / 2 where =^-1+1/2•at,
'\A' ^,1+1/2 • ^i-1/2 2^1+1/2 (A.4) ( Q ?>
0) m < 0) The discrete SaintVenant equations are solved t h r o u g h a t r i diagonal f o r m a t , i f Eq. (A.4) is substituted i n t o (A.3). The t r i d i a g o -nal f o r m a t can t h e n be r e - f o r m e d i n Eq. (A.5). The last equation is i n t e n d e d to create the change o f structure f l o w zl Qc as the control i n p u t f o r o p t i m i z a t i o n . By m u l t i p l y i n g the inverse o f the f i r s t m a -t r i x o n b o -t h sides o f Eq. (A.5), i -t gives -the linear -t i m e v a r y i n g state-space m o d e l f o r m a t ." l l
«li
0 0 0 0 0 As • ' / - l , ( - 2 0 ^1-1,1-11
C2 At 0 0 0 0 0 0 0 At1
rf+1, rk+1 rt<+l C/ 0 b,_2 - ö , _ , Kl (A.5)w h e r e T„,i is the t o p w i d t h at i t h discrete point, Qin is the u p -stream i n f l o w , w h i c h is assumed as k n o w n disturbance.
-At. A X T ; . , A f ,k rk l / " i - l / 2 . ! = 2, 1
= 1,
, 1 - 1 a i =l +
-AxT," < i - i / " , - i / 2 + A | i / u ; ' w , i ' " i - l / 2 (A.1) / u J! = 1,
i+1/2 1 , . . , , ; + 1 / 2 K y 2 -• / "i+\/2 + V\ r+1/2One o f the c o n t r o l goals is to m a i n t a i n the d o w n s t r e a m w a t e r level to the target level (^Q). Therefore, the state is the error bet w e e n bethe bet w o levels w r i bet bet e n as: x = &lbet;; ~ q^. Because o f bethe s y m m e -t r y o f -the f i r s -t m a -t r i x on -the l e f -t side o f Eq. (A.5), q i n -t h e s-ta-te vector can be s i m p l y replaced by x and the e q u a l i t y f u n c t i o n re-mains (all the w a t e r levels along the reach subtract the same target level
qt)-References
111 Wahlin BT. Performance of model predictive control on ASCE test canal 1. J Irrig
Drain Eng 2004;130(3):227-38.
[2] Willems P, Blanco TB, Chiang PK, Cauwenberghs K, Berlamont J, De Moor, B. Evaluation of river flood regulation by means of model predictive control. In: 4th international symposium on flood defence. Managing flood risk, reliability and vulnerability. 2008.
|3] van Overloop PJ. Model predictive control on open water systems. Delft: lOS Press; 2005.
|4] Rodellar J, Gomez M, Bonet L. Control method for on-demand operation of open channel flow. J Irrig Drain Eng 1993;119(2):225-41.
|5] Negenborn RR, van Overloop PJ, Keviczky T, De Schutter B. Distributed model predictive control of irrigation canals. Netw Heterogen Media (NHM) 2009;4(2):359-80.
[6] Barjas Blanco T, Willems P, Chiang PK, Haverbeke N, Berlamont J, De Moor B. Flood regulation using nonlinear model predictive control. Control Eng Pract 2010;18(10);1147-57.
[7] Camacho E, Bordons C. Nonlinear model predictive control: an introductory review. Assessment and future directions of nonlinear model predictive control; 2007. p. 1-16.
[8] Ding Y, Wang SSY. Optimal control of open-channel flow using adjoint sensitivity analysis. J Hydraul Eng 2006:132:1215.
|9] Schwanenberg D, Verhoeven G, Raso L. Nonlinear model predictive of water resources systems in operational flood forecasting. In: 9th international conference on hydroinformatics. HlC 2010, Tianjin, China; 2010.
[10] Xu M, van Overloop PJ, van de Giesen NC. On the study of control effectiveness and computational efficiency of reduced Saint-Venant model in model predictive control. Adv Water Res 2010;34(2):282-90.
[11] Schittkowski K. NLPQ.L: A FORTRAN subroutine solving constrained nonlinear programming problems. Ann Oper Res 1986;5(l):485-500.
[12] Chow VT. Open-channel hydraulics. New York: McGraw-Hill Book Co. Inc; 1959.
[13] Stelling GS, Duinmeijer SPA. A staggered conservative scheme for every Froude number in rapidly varied shallow water flows. Int J Numer Meth Fluids 2003;43:1329-54
[14] Simon D. Optimal state estimation: Kalman, Hoo and nonlinear approaches. New Jersey: John Wiley & Sons, Inc.; 2006.
[15] Gill PE, Murray W, Wright MH. Practical optimization. London, UK: Academic Press; 1981,
]16] Han SP, A globally convergent method for nonlinear programming. J Optim Theory Appl 1977:22:297.
[17] Powell MJD. The convergence of variable metric methods for nonlinearly constrained oprimization calcularions. Nonlinear Program 1978;3:27-63, [18] Powell MJD, A fast algorithm for nonlinearly constrained optimization
calculations, Numer Anal 1978;630:144-57,
]19] Bonnans JF, Gilbert JC, Lemarechal C, Sagasrizabal CA, Numerical opUmizadon: theoretical and practical aspects, Springer-Veriag, New York; 2006.
[20] Sun W, Yuan Y. Optimization theory and methods: nonlinear programming, vol, 1, Springer; 2006,
[211 Clemmens AJ, Kacerek TF, Grawitz B, Schuurmans W. Test cases for canal control algorithms. J Irrig Drain Eng 1998;124(l):23-30.
