- Email: [email protected]
A linear programming regulator applied to hydroelectric reservoir level control
Automatica, Vol. 22, No. 5, pp. 533 541, 1986. Printed in Great Britain. i'
0005 1098/86 $3.00 + 0,00 Pergamon Journals Ltd. 1986 International Federation of Automatic Control
A Linear Programming Regulator Applied to Hydroelectric Reservoir Level Control* PER-OLOF GUTMANt
For linear systems with linear state and control constraints, the regulator problem can be formulated as a linear programming problem. This principle is illustrated by successful applications to a laboratory process and to the reservoir level control of a hydroelectric power station. Key Words--Closed loop systems; computer control; control applications; control non-linearities; feedback control; level control; linear programming; linear systems; non-linear control systems; optimal control; power station control.
Abstract--For linear dynamical systems with linear state and
using only the first p o r t i o n of the c o m p u t e d c o n t r o l trajectory, a feedback c o n t r o l is achieved, the O p e n L o o p O p t i m a l F e e d b a c k ( O L O F , Dreyfus, 1964). V a r i o u s aspects of the L P - O L O F idea, especially for the m i n i m u m time p r o b l e m have been studied over the years ( C a n o n et al., 1970; Bashein, 1971; C a d z o w , 1974; R a s m y a n d H a m z a , 1975; Kolev, 1978; C h a n g a n d Seborg, 1981; Vlieger et al., 1981, 1982). In water resource systems, L P has been used for various design, p l a n n i n g a n d open l o o p feedback c o n t r o l purposes, m a i n l y in the long a n d m e d i u m time range, e.g. M a a s s (1962), Revelle et al., (1968, 1969), L o u c k s (1968), N a y a k a n d A r o r a (1971). O t h e r o p t i m a l m e t h o d s , like d y n a m i c p r o g r a m m i n g have been suggested for short term o p e r a t i o n also, e.g. L a r s o n a n d K e c k l e r (1969), where a 24-h o p e r a t i o n with l - h s a m p l i n g interval is suggested. In S a c h d e v a (1982) an extensive bibliog r a p h y on these a n d related topics is presented. It does n o t seem, however, that an L P r e g u l a t o r o p e r a t i n g in the o p e n l o o p feedback m o d e has been used in a h y d r o e l e c t r i c water reservoir system in the short time range in place of a c o n v e n t i o n a l controller (e.g. P I D controller) for the i m m e d i a t e c o n t r o l of some variable. F o r a chemical p l a n t with constraints on the c o n t r o l inputs, K n u d s e n (1975) i m p l e m e n t e d an L P - O L O F regulator. In this p a p e r (condensed from G u t m a n , 1982a, b, 1983a, b) an L P - O L O F r e g u l a t o r for linear systems with b o t h c o n t r o l a n d state constraints is i m p l e m e n t e d to control, on line, a d o u b l e water t a n k l a b o r a t o r y process a n d the water level of a hydroelectric p o w e r station reservoir. The cont r i b u t i o n of this p a p e r is believed to be the testing of the L P - O L O F r e g u l a t o r to c o n t r o l real life processes with c o n t r o l and state constraints.
control constraints the regulator problem can be formulated as a linear programming problem. A regulator, built around a standard LP program and operating in the Open Loop Optimal Feedback (OLOF) fashion, is presented. The LP-OLOF regulator was implemented on a VAX 11/780 computer to control, in real time, a double water tank laboratory process, and the water level of a hydroelectric power station reservoir. The reservoir control experiment showed that even with an assumed simple process model, satisfactory performance was achieved. In particular it was beneficial that the LP-OLOF regulator allows dynamic changes of the water flow bounds corresponding to the number of generators available at different operating conditions, and that the regulator can predict the water level. 1. INTRODUCTION LINEAR PROGRAMMING as a tool for a u t o m a t i c c o n t r o l has e v o k e d interest for at least 25 years. E a r l y c o n t r i b u t i o n s include M a n n e (1960), who used linear p r o g r a m m i n g (LP) to c o m p u t e state feedback c o n t r o l laws, a n d M a a s s et al. (1)62) where L P was investigated as a m e t h o d to determine design a n d seasonal o p e r a t i o n of m u l t i - p u r pose m u l t i - r e s e r v o i r water resource systems. Z a d e h (1962) a n d P r o p o i (1963) noticed that the m i n i m u m time p r o b l e m for linear discrete systems with control c o n s t r a i n t s can be f o r m u l a t e d as a series of L P problems. P r o p o i (1963) p r o p o s e d that the L P p r o b l e m s be solved in each s a m p l i n g interval. By * Received 10 March 1984; revised 7 October 1985; revised 7 February 1986. A preliminary version of this paper was presented at the 3rd IFAC/IFIP Symposium on Software for Computer Control which was held in Madrid, Spain during October 1982. The published proceedings of this IFAC meeting may be ordered from Pergamon Books Limited, Headington Hill Hall, Oxford, OX30BW, England. This paper was recommended for publication in revised form by Associate Editor Joe H. Chow under the direction of Editor H. Austin Spang III. 1"Israel Electro-Optics Industries, P.O. Box 1165, Rehovot 76110, Israel. 533
534
P.-O. GUTMAN
The paper is organized as follows. Section 2 briefly describes the considered control problem and various methods to solve it, including the LPO L O F approach. In Section 3 the control of the laboratory process is found and in Section 4 the hydroelectric water reservoir test is reported. Sections 5 7 cover the discussion, acknowledgement and references, respectively.
2. PRELIMINARIES Assume that a plant is described by a linear discrete model, with constant control and state constraints (t + 1) = ~x(t) + Fu(t) + f ( t )
(2.1a)
au(j) <~ uj4t) <~ [~,(j),j = 1,2 . . . . . m;allt >~ t o
(2.1b)
ax(j) <~ xi(t) <<.[3x(i),i = 1,2 . . . . . n;allt >~ to,
(2.1c)
where f is a measured, b o u n d e d disturbance, x is an n-vector, u an m-vector, and ~ , F matrices of appropriate dimensions, c~and/~ define the controland state constraint sets. It is assumed that the whole state is measured (or estimated) at each sampling instant. More general linear time invariant constraints can also be introduced. The objective is to design a regulator that takes the state (or output) from its initial value to a prescribed stationary point, which will define the origin in the state space. It is assumed that the state and control constraints are admissible, i.e. the initial state is such that it is possible to stabilize (2.1) to the origin without violating the constraints. In G u t m a n (1984) necessary and sufficient conditions for admissibility are given when f = 0, and a variable structure linear state feedback controller is suggested. A straightforward way to solve the above control problem is to design a linear controller that never violates the constraints. Linear quadratic design may be used, or model following, i.e. the set point is changed in a modest way (e.g. a ramp reference input). Simulation must be used to check that these methods work. The c o m m o n method in practice is to design a linear controller disregarding the constraints and saturating the control input. The region of stability may be very small. For open loop stable plants Astr6m (1971) proposes a stabilizing relay, and Shapiro (1972) determines the matrix L so that u(t) = sat(Lx(t)) is stabilizing. Under certain conditions there exists a stabilizing u(t)= sat(Lx(t)) also for open loop unstable plants (Gutman, 1982a). Glattfelder (1974) and others have proposed the method of m a x - m i n selectors, whereby a set of controllers is used. W h e n a state constraint is violated, the control is switched to a controller
whose aim is to remove the state from the boundary. Although there is no proof of stability (the states may oscillate from one b o u n d a r y to another) the m e t h o d has become quite popular in e.g. the power industry. An elegant way to design a controller for (2.1) is to solve an optimal control problem with an optimization criterion reflecting the desire to stabilize the plant. In general, however, one gets an open loop control solution from one initial state, not the desired feedback control. To construct a feedback control, one may solve the optimal control problem for a great n u m b e r of initial conditions and then set up a look-up table. Alternatively, solve the optimal control problem on line at every sampling instant, and deliver the c o m m a n d signal pertaining to the current sampling interval.* This is the Open L o o p Optimal Feedback scheme (OLOF). Since the plant and constraints in (2.1) are linear, it is possible to use L P to solve the optimal control problem. However L P does not allow a free endtime. A free end-time problem can be solved by solving a sequence of fixed end-time problems, iterating over the end-time. Here the end-time is called "the time horizon r" and the end-time that solves the free end-time problem is called "the optimal time horizon Z'opt". Since it is assumed above that there exists a stabilizing control sequence, there exists at least one r for which the L P problem has a solution. So, at each sampling instant two or more optimization problems are solved for consecutive values of the time horizon z to find the predicted u(0), u(l) . . . . . u(z - 1), and x(1), x(2) . . . . . x(r) where r = 0 is present time. The optimization criterion is min ~. ( r l m a x [ - x i ( t )
+ 7.~(i),0]
i , t <<. z
+ wimax[xi(t ) - fix(i),0] + ~ cilxi('c)t,
(2.2)
i
where the vectors r and w are b o u n d a r y region penalty weights, ~'x and fix define the b o u n d a r y regions, and c is the terminal penalty weight. The optimal time horizon ropt is assumed to be found when Ix(r)la ~ ~', diixi(r)l ~< *:.
(2.3)
i=1
When solving the free end-time problem, it is essential that the iterations continue until topt is reached. A shorter final r might be harmful because stabilization to the origin might never be achieved. A longer final r can also be harmful in the O L O F setting, since the (non-unique) control obtained for the current sampling interval might counteract the control objective. However, during real time *The estimated computation time can be introduced as a constant delay in (2.1a).
A linear programming regulator operation it might happen occasionally that one sampling interval is not sufficient to solve the number of LP problems needed to find Topt. In such cases, the control suggested by the solution for the latest computed z is used. This is called "the suboptimal control". A standard linear programming package, PRIMAL, supplied by the Institute of Applied Mathematics, Stockholm (Holmberg, 1981) was used. The base inverse matrix is stored in product form, and reinversion is done when the multiplication becomes too time consuming, see Orchard-Hays (1968). The L P - O L O F controller was implemented on a VAX 11/780 computer, and the simulations done with the simulation language S I M N O N (Elmqvist, 1975). More details about the control algorithm are found in Gutman (1982a, b). 3. C O N T R O L
OF A LABORATORY
PROCESS
The LP regulator was used to control a laboratory process consisting of two water tanks with free outlets, as drawn in Fig. 1. The water heights were measured. The control objective was to get h2 to a reference value, h2ref , a s fast as possible, without causing overflow in either tank. The input voltage to the pump was limited. The non-linear pump characteristic was found to be so harmful that a local, analogue PID-feedback around the pump was introduced, whose dynamics were sufficiently fast to be neglected. Figure 2 depicts the block diagram of the process. The double tank system is non-linear. The time constant (from u to h2), when filling is about 125 s;
but only about 60 s when emptying. The sampling interval was set to 12s, because it is a suitable fraction of the shorter time constant, and because 12 s suffice as computation time for most reference changes. Dimensionless units were introduced: O~u~l,
~ Tank I
l
tocho signal I
= -A input voltagg
~
I
Pump [
]
Tonk
water FIG. I. The double tank process: physical set-up,
I
b2ref
FIG. 2. The double tank process; block diagram.
O~hl~l
,
O~h2~l,
(3.1)
u = 0 means no water flow, u = 1 means maximum flow. For ht and h 2 , 0 means that a tank is empty, 1 that it is full. The following continuous and sampled models were obtained from some step input experiments.
h(t) = \(--001970.0178 --0.01290) h(t) + _t/0"0263~_ LO )u(t) (3.2) "~ h(t + 12) = \0.175(0'7900.8570) h(t) + /'0.282 L0.0296/u(t) (3.3) where h(t) = [ht(t)h2(t)] v. This model was used in the LP regulator, with ULp ~ U - - Uref,
X 1 =- h l r e f ,
(3.4)
X 2 = h 2 -- h2ref
with hlref and u,,f computed from the stationary solution of (3.2), given h2ref. The input and state constraints are given in (3.1). In (2.2) and (2.3) e=0.05,
ftow
535
7x = [ - - h l r e f
r=w=[11], 0 ] T,
c=d=[l
6 x = [-hlref
Jr-
4],
0 . 9 0 ] x, (3.5)
were chosen. Since h2ref = 0 is an acceptable reference, which means ht = 0 in stationarity, there was no boundary penalty for low tank 1 levels. Therefore, y:,(1)= -hlref. The values of 7x(2) and 6x(2) mean that Ih2 - h2refl is included in the loss for all t, not only at the time horizon z. This is not the intended way to define a boundary region but it was decided on after a few preliminary experiments. The behaviour for small deviations from h2ref w a s especially improved. Moreover, when the sampling interval is not enough to find the optimal time horizon, it is good to penalize the deviation from h2ref along the trajectory to ensure that the control leads the system in the right direction. Since the second tank is essentially a first order system, this type of penalty function will not affect the time optimality of the LP solution. With 7x(2) 4:0 and 6,(2) :~ 0, the control tended to become a dead-beat control for small deviations from h2ref. This caused oscillations since the model
536
P.-O. GUTMAN
0
0
t000
0
1000
~___-f
.~__f
I reference
I/
i\ lo'oo
0
loICOmputing time [sees} in ir eoch sompling interval
O,L 0
,
10~)0
FIG. 3. The double tank process: u(t), hl(t ), hzref(t) ( = reference) and h2(t) during a real life experiment and the computation time in the simulated experiment.
was not perfect. It was also found that a minimum time horizon of two sampling intervals was beneficial for the local control, in addition to the above boundary region. The tanks were initially empty during the real life experiment, h2 was to follow h2ref , see Fig. 3. No outer disturbances affected the system. The computation times are perceptible in the figure. At reference value changes they were up to 4 s, except at the last reference value change, when the computations were interrupted after 10 s, and the suboptimal solution, u = 0, corresponding to • = 23, was used. The interruption was made because it was estimated that the next iteration would make the total computation time longer than 12 s. In "steady state", the computation times were less than 1.5 s. At reference value changes, the time horizon was restarted at 1. For reference value changes upwards, r roughly reflected the settling time, while it was too long for downward changes. This explains the overshoot at downward reference changes. The computation times at the reference changes were: for ~ = 4 , 0.23s; for • = 9 , 1.1s; for ~ = 11, 1.8 s; for ~ = 15, 3.4 s; and for r -- 23 (suboptimal), 10 s. These data suggest that the computation time is approximately proportional to r2 ln~. Notice that the state constraints are not violated, and that the local control is smooth. As a comparison, this LP regulator was simulated against the linear model (3.2). The computation times were included in the
simulation, but not in the sampled model. See Fig. 4. Notice the computation times in Fig. 3, which are roughly equal to the real life computation times. There are two main deviations between real life and simulation: in real life there is a slight stationary error in h2; the real life downwards dynamics are faster than the simulated. These discrepancies are due to the use of the simple linear model (3.3). The conclusions from the laboratory tank experiment, and supporting simulations, were, briefly: the L P - O L O F regulator can be used even if the computation time occasionally is longer than the sampling interval; the regulator functions well even if an approximate linear model of the non-linear
simulofed
I
/
/
h2(f)
//
0
0
1000
FIG. 4. The double tank process: h2(t) in the real life and simulated experiments.
A linear programming regulator plant is incorporated into it; the performance is sensitive to the tuning of the penalty function weights; the state constraints were rarely violated and the control constraints were never violated. The conclusions encouraged the author to find a real life process with suitable dynamics and essential control and state constraints on which the L P - O L O F regulator could be tested. The level control problem of the reservoir of the Bergeforsen hydroelectric station was thought to be suitable. The laboratory tank process could be seen as a very simple model of the reservoir. Therefore the approach to solve the Bergeforsen problem was similar: identification and modelling of a linear model, simulation, tuning of the penalty function weights and experimentation. 4. T H E B E R G E F O R S E N LEVEL C O N T R O L EXPERIMENT
The Bergeforsen reservoir, consisting of the lowest 55 km of Indals/ilven near Sundsvall, Sweden has at its upper end the Jiirkvissle power station and at its lower end the Bergeforsen power station. The flow at Bergeforsen ranges from ~ 5 0 m 3 s-1 to ~ 840m 3 s - i s and almost all the water enters from J/irkvissle. The Water Court decision, governing the operation of the power plants, states that the water level at K/ivsta, 14.3 km upstream from Bergeforsen "must coincide as close as possible but not exceed 23 m above sea level". See Fig. 5. At present the J/irkvissle power station is run according to the grid power demand while the four turbines at Bergeforsen are controlled manually to keep the K/ivsta level correct. Typically the J/irkvissle flow will be at a high level in the day and at a low level at night. The operators at Bergeforsen get oral information over the phone in advance of major changes in the J/irkvissle flow. The three upper curves of Fig. 7 represent normal operation over 75 h. A block diagram model of the automatic system is suggested in Fig. 6, with the J~irkvissle flow u 2 as a measured disturbance (feedforward), and the Kfivsta level y as the output feedback signal. The Bergeforsen flow Ul is the control signal. The control objective is to keep Jiirkvissle
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Kavsto level
Bergeforsen r---~----~ flow u~ ~ I Controller ~ Delay t - - ' - - ' ~ Dynamics I ~.___~r___a (m~/s) i , i i
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
y(mS/s)
.
.
.
.
.
.
.
.
J
FIG. 6. Block diagram of the Bergeforsen reservoir level control system.
y ~< 23 m, but as close as possible to 23 m in order to concur with the Water Court decision, and to utilize the reservoir capacity maximally. A transport delay is assumed from the inputs ul, u 2 followed by linear dynamics. Such a model is approximate since the reservoir dynamics are best described by non-linear partial differential equations. Previous data and physical insight suggest that the delay from ul should be 20min, the delay from u2 c a 1 h, and that the dynamics should include an integrator and some wave phenomenon. An identification experiment was performed, where especially u~ was varied randomly around the operating level, see Fig. 7. The sampling interval was set to 5 rain and necessary anti-aliasing prefiltering was done with analogue and digital filters. Correlation analysis, least squares and maximum likelihood methods of the interactive program Idpac (Wieslander, 1980) were used to identify the delays and dynamics. The identified model output is seen in Fig. 7. This model was modified when the L P - O L O F regulator was simulated against a computerized PDE model ot the reservoir, Routes, developed at the Swedish State Power Board. The final model was in state space form: x ( t + 1) = ~ x ( t ) + FUl(t ) + F(t)
(4.1a)
with x(t) = (y(t)ul(t ul(t -
- . 1)ul(t - 2)
3)ut(t - 4)ul(t - 5)
(4.1b)
u l ( t - 6)Aul(t - 2) Aul(t - 2)) r, zXul(t) = u l ( t + 1) -- u~(t).
,e..,[m]
< <
~ m31s Jarkvissle flow
.
Kavsta
[m3,s]
,,o.
537
5 5 Km
FIG. 5. The Bergeforscn reservoir.
Bercjeforsen
) )
(4.1c)
538
P.-O. GUTMAN
o
zoo
600
Somple 800
.oo
,oo
,oo
,oo
600
800
400
no,
,ooo
oJ[m]
ooT-- . . . . . . . . . . . . . . .
500 ; f ~ _ _ _ _ . v S ~
, Som~ no. [
o
,ooo
LK[ml
zz.91.-
_ '...
o
. .
Sample no.
~
20°
,o®
YM 0
.
0
.
200
.
.
400
600
l
800
I000
FIG. 7. Identification experiment, measurements for Bergeforsen flow (QB), J/irkvissle flow (QJ), KS.vsta level (LKI and idcntified model output for K/ivsta level with mean and linear trend subtracted (f'). 4.7451'10 -5
0
0
0-
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 0
0
1
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0 0
0 0
0 0
0 0
0 0
0 I
0 0 1
1
0
0
--5.7106" 10 -5
0
0
0
0
0
1
0
0
0
1
0
0
0
0
0
0
0
-1
0 0
-1
F(t)
= ((1.02-lO-Suz(t - 9)
1
-[-.f(uz,t
- 0 . 0 5 4 5 " 10 -5
i n p u t experiment) b/2o = 3 7 0 m 3 s - 1.
1)) 0 0 0 0 0 0 0 0)I",
+
(4.1e) f(u
2, t) =
lO-5(6.49(u2(t -
tl2(t -
15)) - 5.71(Uz(t - 15) - U z ( t -
- 0.78(uz(t - 16) --
3.8(U2(t
--
--
b~2(t
--
16))
(4.1f)
With the integrator emphasized, the model (4.1) is written y ( t + 1) = y ( t ) - 1.02" 10-5(ul(t -- 3)
17))
20) - u 2 ( t - 21))
+ 1.1(u2(t - 21) - u z ( t - 22)) + 2.7(Uz(t - 22) - u 2 ( t - 23))),
where u l , u2, y stand for Bergeforsen flow c o m m a n d input, JS.rkvissle flow a n d K&ivsta level a r o u n d their respective, n o m i n a l values. Ulo = 3 9 0 m 3 s - l Yo = 22.95 m (changed on line to 22.94 in step
(4.2)
The difference between Ulo a n d Uzo represents the average effect of side-flows into the reservoir. Aul was i n t r o d u c e d as two extra states in order to allow piecewise linear penalties on flow changes (see Fig. 8d). F ( t ) represents the feedforward term from u2, and f ( u 2 , t + 1) the wave generated by u2. Note that the s t a t i o n a r y value of f is 0.
14)
-
(4.1 d)
-
u2(t
-
9))
+ (terms representing waves) (see (4.1e) a n d note that
A linear p r o g r a m m i n g r e g u l a t o r
539
5
~b(l,i) = - 1 . 0 2 . 1 0 - 5 in (4.1d)).
The c o n s t r a i n t s a n d penalties are illustrated in Figs 8(a)-(d). Notice that high water levels are more penalized t h a n low (cf. W a t e r C o u r t decision), that the flow b o u n d s were c h a n g e d on line d e p e n d i n g on the n u m b e r of g e n e r a t o r s and that the introd u c t i o n of two Au t states m a d e it possible to shape the flow change p e n a l t y function, so that small changes are cheap, large changes costly a n d very large changes f o r b i d d e n (since they have d i s a s t r o u s effects on the delta d o w n s t r e a m from Bergeforsen). The m i n i m u m time h o r i z o n was set to 10, so that the current u 2 m e a s u r e m e n t w o u l d always enter (4.1), the m a x i m u m time horizon was set to 23, and the m a x i m u m C P U time before a s u b o p t i m a l c o n t r o l was used, 200s, The final c o n t r o l experiment e n c o m p a s s i n g 190 samples is shown in Fig. 9. T h e n u m b e r of g e n e r a t o r s was decided a c c o r d i n g to the p l a n n e d J/irkvissle flow, and the u~ b o u n d s were set accordingly. At 03.03 hours the Kfivstalevel reference Yo was c h a n g e d to 22.94 m. It is clear that y was c o n t r o l l e d satisfactorily, in general on the low side of Y0, with an acceptable oscillative a m p l i t u d e (cf. Fig. 7 when the control was manual).
i=3
T o tune the p a r a m e t e r s of the L P - O L O F controller, it was f o u n d to be beneficial to simulate it against m o d e l (4.1), i.e. the m o d e l internal to the controller. This is called i n t e r n a l s i m u l a t i o n . M o r e over this m o d e l can be used to p r e d i c t future o u t p u t s d u r i n g o p e r a t i o n . W h e n the p r e d i c t i o n time o r the time h o r i z o n z of the c o n t r o l a l g o r i t h m is longer t h a n 10 s a m p l i n g periods, then future, as yet u n m e a s u r e d u 2 values enter (4.1). u2 was then predicted in the simplest possible way: U2('C + g) = Uz(t), Z > 0, where u2(t ) = latest m e a s u r e m e n t of u 2. In a refined controller, future k n o w n changes of u 2 can be a c c o u n t e d for. T h e final p a r a m e t e r set was as shown in Table 1, a n d % = Ctx(i ),
i = 1,2 . . . . . 7
ft, = fix(i),
i = 1,2 . . . . . 7.
(4.3)
O n the first occasion, when 3 turbines were used, the u p p e r u l - b o u n d was 230, the second time it was 240. It can be shown that for limited u2 a n d Au 2, (4.2) a n d T a b l e 1 are admissible c o n s t r a i n t sets. In (2.2), (2.3), ci
=
10;
d i = ci;
rt
=
200;
to i = 1000;
c i = 0,
i=
i = 2 . . . . . 9;
1. . . . . 9;
Ti = 0 ,
i = 2 . . . . . 7;
O,)i ~ ri~
i = 2 . . . . . 9;
r 8 = 0.1,
r9
=
0.01;
(4.4)
= 0.05. TABLE 1.
S t a t e no. (5.1)
~tx(i) (2.1c)
)'x(i) (4.2)
6~(i)
flxll)
(4.2)
(2.1c)
Level Flow
1 2..... 7
zXu* Au,*
8 9
-0.15 - 320 - 290 -100
-0.004 - 320 - 290 -15 0
0.004 20 230/240 15 0
0.15 20 230/240 100 100
Eqn. no:
* nominal
value
-100 =
Operating condition: no. of turbines 2 3
0
Intermediate
Terminal level penalty
level
penalty
10o
I
' I
-0.15
stopping criterion
[m] 0.15
toleronee
FI(;. 8(a), P e n a l t y f u n c t i o n c o m p o n e n t s : t e r m i n a l level penalty.
I
-0.15
I
0
QI5
FIG. 8(b). Penalty function components: intermediate level penalty.
;-
540
P.-O. GUTMAN One exception is the overshoot at 08.45 of 3.4cm which did not fall within the unofficial tolerance at 3 cm. It is believed that the reason is the insufficiency of the model; in particular, the time delay from u~ should have been 35 min instead of 20 min. The mismodelling is probably due to weak identifiability at the delay. In the identification experiment (Fig. 7) u~ changed only a couple of times. In the judgement of the operators who run the plant manually, the L P - O L O F control was satisfactory, corresponding to an average operator. The advantages were seen as (1) the possibility to change Bergeforsen flow bounds on line; (2) the flow bounds were not violated; (3) the prediction of future K/ivsta levels and control inputs, the latter making it possible to govern the flow (suboptimally) for a few sampling intervals even when the computer went down (this happened at one occasion); (4) the feedforward from the J/irkvissle flow; (5) the nonsymmetrical penalty on y and the possibility to restrict flow changes. All in all, with somewhat more exact modelling the L P - O L O F regulator seems viable for the reservoir level control problem.
Flow penalty '(2 turbine case)
I -300
I -200
I -I00
I I00
)
FI(,. 8(c). Penalty function components: flow penalty. Flow c h a n g e 'penalty
l
,o
5. C O N C L U S I O N S
Of course, few general conclusions about the LPO L a F regulator can be drawn from two experiments. The most obvious advantage over most alternative methods is that the state and control constraints are taken into account when the control is computed, and that a criterion reflecting the design goal is optimized. Some other general points arc discussed in G u t m a n 1982, 1983). Here some specilic remarks are made concerning the on-line
3/S/5rain] I
-Io0
o
)
I0o
I"1(;. 8(d). Penally, function components: flow change penalty.
,[m3/s] f l o w
r
600 500. . . . "~"
QB-bounds
400-
J
. . . . . . . . .
....... ........
500200I00time[hours]
t 02.15
I
2015 ,[m]
2215
0015
I
0415
I 0915
I
0615
I
level
2297
:l :' "'-..
22 95
:
,'.. ", ......: ......L K
,-.Y"., \,.
\
22 95-
,/.Y-
_
:""... / ...... • :
: ."
I
:, L K r e f ., -"...................
..
. :,.:
""":; 2015
)
1015
2215
O015
L,
'." A .........-
............
time[hours] 02.15
0415
06.15
0815
10.15
I't~;. t.L l)iagran~ of Imal cxperimcm: Bergcforsen I1~',', ~QB), Jiirk,,isslc flow (QJ). Kfivsta level H_KI and Kii,,sla rclk:rcnce level (kKrcf).
A linear p r o g r a m m i n g r e g u l a t o r e x p e r i m e n t s . T h e r e g u l a t o r p r o v e d itself r o b u s t f r o m t h r e e p o i n t s o f view. (1) I n spite o f the a s s u m e d h i g h l y a p p r o x i m a t e d l i n e a r m o d e l s of the n o n l i n e a r plants, it p e r f o r m e d well. (2) It w i t h s t o o d well o c c a s i o n s w h e n the c o m p u t a t i o n t i m e s were l o n g in c o m p a r i s o n to the s a m p l i n g interval, by u s i n g s u b o p t i m a l c o n t r o l s ; a n d a n o c c a s i o n of computer shutdown when previously computed p r e d i c t e d c o n t r o l s w e r e used. (3) It is sensitive to t h e t u n i n g of the p e n a l t y f u n c t i o n weights. T h e B e r g e f o r s e n e x p e r i m e n t s h o w e d t h a t o n e of the m a i n a d v a n t a g e s of the L P - O L O F r e g u l a t o r is the p o s s i b i l i t y o f h a v i n g h a r d b o u n d s o n the c o n t r o l variables and their derivatives, which are taken into a c c o u n t w h e n the c o n t r o l a c t i o n is c o m p u t e d (cf. w i t h s a t u r a t o r s a c t i n g o n the a l r e a d y c o m p u t e d control). O t h e r advantages include dynamic c h a n g e s of b o u n d s , s k e w e d p e n a l t y functions, pred i c t i o n s o f f u t u r e p r o c e s s values a n d the possibility of t u n i n g the L P - O L O F by s i m u l a t i n g it a g a i n s t itself. Acknowledgements--The Bergeforsen experiment was performed
under the aegis of Kockumation AB, MalmS, (thanks to J. Tuszynski) and the Swedish State Power Board (thanks to A. Bruse and G. Fabricius). The report was completed during the author's stay at the Technion-Israel Institute of Technology, which was made possible by generous grants from the C. F. LundstrSm Foundation, the Swedish Board for Technical Development, the Helge Ax:son Johnson Foundation and the Per Westling Foundation of Lund University. The author thanks the reviewers, in particular Reviewer # 2, who suggested several relevant references on water systems control.
7. REFERENCES AstrSm, K. J. (1971). Olinj~ira system (in Swedish). Lund Inst. Techn., TLTH/VBV, Lund, Sweden. Bashein, G. (1971). A simplex algorithm for on-line computation of time optimal controls, lEE Trans. Aut. Control, AC-16, 479-482. Cadzow, J. A. (1974). Minimum amplitude control of linear discrete systems. Int. J. Control, 19, 765. Canon, M. D., C. D. Cullum and E. Polak (1970). Theory of Optimal Control and Mathematical Programming. McGrawHill, New York. Chang, T. S. and D. E. Seborg (1983). A linear programming approach for multivariable feedback control with inequality constraints. Int. J. Control, 37, 583-597. Dreyfus, S. (1964). Some types of optimal control of stochastic systems. S l A M J. Control, 2, 120. Elmqvist, H. (1975). Simnon An interactive simulation program for nonlinear systems. Rep. 7502, Dept. Aut. Control, Lund Univ., Sweden. Glattfelder, A. H. (1974). Regelungsysteme mit Bergrenzungen. Oldenburg, Mi.inchen. Gutman, P.-O. (1982a). Controllers for bilinear and constrained linear systems. Ph.D. dissertation, Dept. Aut. Control, Lund Univ., Sweden, CODEN:LUTFD2/TFRT- 1022)/1-133/ (1982).
AUTO 22:5-C
541
Gutman, P.-O. (1982b). On-line use of a linear programming controller. 3rd IFAC/IFIP Symp. Software Comp. Control, Madrid, Spain, October 5-8. Gutman, P.-O. (1983a). On the linear programming OLOFregulator. Technion-Israel Inst. Tech., Haifa, Israel, EE Pub. No. 469. Gutmam P.-O. (1983b). Test of the OLOF-regulator for water reservoir level control. Technion-Israel Inst. Tech., Haifa, Israel, EE Pub. No. 468. Gutman. P.-O. (1984). Admissible sets and feedback control for discrete time linear dynamical systems bounded controls and states. Proc. 23rd IEEE Conf. Decis. Control, December 1214, Las Vegas, Vol. 3, pp. 1727-1731; IEEE l~'ans. Aut. Control, AC-31, pp. 373-376 (1986). Holmberg, G. (1981). PRIM2-ett subrutinpaket som 15ser LPproblem (in Swedish). UTM-Prog. nr. 50. Uppdragsgruppen f. tillfimpad matematik. Royal Inst. of Techn., Stockholm, Sweden. Kolev, L. V. (1978). Minimum-fuel minimum-time control of linear discrete systems. Int. J. Control, 27, 21-29. Knudsen, J. K. H. (1975). Time-optimal computer control of a pilot plant evaporator. Proc. IFAC 6th World Congress, Part 2, pp. 45.3/1-6. Larson, R. E. and W. G. Keckler (1969). Application of dynamic programming to the control of water resource systems. Automatica, 5, 15-26. Loucks, D. P. (1968). Computer models for reservoir regulation. J. Sanit. Eng. Div., Amer. Soc. Civil Eng., 94(SA-4), 657-669. M aas, A. et al. (1962). Design of water resource systems. Harvard University Press, Cambridge, Massachusetts. Manne, A. S. (1960). Linear programming and sequential decisions. Manag. Sci., 6, 259-267. Nayak, S. C. and S. R. Arora (1971). Optimal capacities for a multi-reservoir system using the linear decision rule. Water Resource Res., 7, 485 498. Orchard-Hays, W. (1968). Advanced Linear Programming Computing Techniques. McGraw-Hilk New YorL Propoi, A. I. (1963). Use of L P methods for synthesizing sampleddata automatic systems. Aut. Remote Control, 24, 837. Rasmy, M. E. and M. H. K. Hamza. Minimum-effort timeoptimal control of linear discrete systems. Int. J. Control, 21, 293-304. Revelle, C., D. P. Loucks and W. R. Lynn (1968). Linear programming applied to water quality management. Water Resource Res., 4, 1-10. Revelle, C., E. Joeres and W. Kirby (1969). The linear decision rule in reservoir management and design I. Development of the stochastic model. Water Resource Res., 5, 767-777. Sachdeva, S. A. (1982). Bibliography on optimal reservoir drawdown for the hydroelectric-thermal power system operation. IEEE Trans. Power Apparat. Syst., PAS-101, 14871495. Shapiro, E. Y. (1972). Suboptimal design of minimal-time feedback regulator for linear time-invariant plants. Sc.D. Dissertation, Dept. Elec. Eng. Comp. Sci., Columbia Univ., New York. Vlieger, J. H. de, H. B. Verbruggen, P. M. Bruijn (1981). A time-optimal control algorithm for on-line computer control allowing bounding control and state variables. Control Comp., 9, 1-6. Vlieger, J. H. de, H. B. Verbruggen, P. M. Bruijn (1982). A timeOptimal control algorithm for digital computer control. Automatica, 18, 239-244. Wieslander, J. (1980). Idpac commands--user's guide. Dept. Aut. Control, Lund, Sweden, CODEN: LUTFDZ/(TFRT3157)/1-108/(1980). Zadeh, L. A. and B. H. Whalen (1962). On optimal control and linear programming. IRE Trans. Aut. Control, 7, 45.
0005 1098/86 $3.00 + 0,00 Pergamon Journals Ltd. 1986 International Federation of Automatic Control
A Linear Programming Regulator Applied to Hydroelectric Reservoir Level Control* PER-OLOF GUTMANt
For linear systems with linear state and control constraints, the regulator problem can be formulated as a linear programming problem. This principle is illustrated by successful applications to a laboratory process and to the reservoir level control of a hydroelectric power station. Key Words--Closed loop systems; computer control; control applications; control non-linearities; feedback control; level control; linear programming; linear systems; non-linear control systems; optimal control; power station control.
Abstract--For linear dynamical systems with linear state and
using only the first p o r t i o n of the c o m p u t e d c o n t r o l trajectory, a feedback c o n t r o l is achieved, the O p e n L o o p O p t i m a l F e e d b a c k ( O L O F , Dreyfus, 1964). V a r i o u s aspects of the L P - O L O F idea, especially for the m i n i m u m time p r o b l e m have been studied over the years ( C a n o n et al., 1970; Bashein, 1971; C a d z o w , 1974; R a s m y a n d H a m z a , 1975; Kolev, 1978; C h a n g a n d Seborg, 1981; Vlieger et al., 1981, 1982). In water resource systems, L P has been used for various design, p l a n n i n g a n d open l o o p feedback c o n t r o l purposes, m a i n l y in the long a n d m e d i u m time range, e.g. M a a s s (1962), Revelle et al., (1968, 1969), L o u c k s (1968), N a y a k a n d A r o r a (1971). O t h e r o p t i m a l m e t h o d s , like d y n a m i c p r o g r a m m i n g have been suggested for short term o p e r a t i o n also, e.g. L a r s o n a n d K e c k l e r (1969), where a 24-h o p e r a t i o n with l - h s a m p l i n g interval is suggested. In S a c h d e v a (1982) an extensive bibliog r a p h y on these a n d related topics is presented. It does n o t seem, however, that an L P r e g u l a t o r o p e r a t i n g in the o p e n l o o p feedback m o d e has been used in a h y d r o e l e c t r i c water reservoir system in the short time range in place of a c o n v e n t i o n a l controller (e.g. P I D controller) for the i m m e d i a t e c o n t r o l of some variable. F o r a chemical p l a n t with constraints on the c o n t r o l inputs, K n u d s e n (1975) i m p l e m e n t e d an L P - O L O F regulator. In this p a p e r (condensed from G u t m a n , 1982a, b, 1983a, b) an L P - O L O F r e g u l a t o r for linear systems with b o t h c o n t r o l a n d state constraints is i m p l e m e n t e d to control, on line, a d o u b l e water t a n k l a b o r a t o r y process a n d the water level of a hydroelectric p o w e r station reservoir. The cont r i b u t i o n of this p a p e r is believed to be the testing of the L P - O L O F r e g u l a t o r to c o n t r o l real life processes with c o n t r o l and state constraints.
control constraints the regulator problem can be formulated as a linear programming problem. A regulator, built around a standard LP program and operating in the Open Loop Optimal Feedback (OLOF) fashion, is presented. The LP-OLOF regulator was implemented on a VAX 11/780 computer to control, in real time, a double water tank laboratory process, and the water level of a hydroelectric power station reservoir. The reservoir control experiment showed that even with an assumed simple process model, satisfactory performance was achieved. In particular it was beneficial that the LP-OLOF regulator allows dynamic changes of the water flow bounds corresponding to the number of generators available at different operating conditions, and that the regulator can predict the water level. 1. INTRODUCTION LINEAR PROGRAMMING as a tool for a u t o m a t i c c o n t r o l has e v o k e d interest for at least 25 years. E a r l y c o n t r i b u t i o n s include M a n n e (1960), who used linear p r o g r a m m i n g (LP) to c o m p u t e state feedback c o n t r o l laws, a n d M a a s s et al. (1)62) where L P was investigated as a m e t h o d to determine design a n d seasonal o p e r a t i o n of m u l t i - p u r pose m u l t i - r e s e r v o i r water resource systems. Z a d e h (1962) a n d P r o p o i (1963) noticed that the m i n i m u m time p r o b l e m for linear discrete systems with control c o n s t r a i n t s can be f o r m u l a t e d as a series of L P problems. P r o p o i (1963) p r o p o s e d that the L P p r o b l e m s be solved in each s a m p l i n g interval. By * Received 10 March 1984; revised 7 October 1985; revised 7 February 1986. A preliminary version of this paper was presented at the 3rd IFAC/IFIP Symposium on Software for Computer Control which was held in Madrid, Spain during October 1982. The published proceedings of this IFAC meeting may be ordered from Pergamon Books Limited, Headington Hill Hall, Oxford, OX30BW, England. This paper was recommended for publication in revised form by Associate Editor Joe H. Chow under the direction of Editor H. Austin Spang III. 1"Israel Electro-Optics Industries, P.O. Box 1165, Rehovot 76110, Israel. 533
534
P.-O. GUTMAN
The paper is organized as follows. Section 2 briefly describes the considered control problem and various methods to solve it, including the LPO L O F approach. In Section 3 the control of the laboratory process is found and in Section 4 the hydroelectric water reservoir test is reported. Sections 5 7 cover the discussion, acknowledgement and references, respectively.
2. PRELIMINARIES Assume that a plant is described by a linear discrete model, with constant control and state constraints (t + 1) = ~x(t) + Fu(t) + f ( t )
(2.1a)
au(j) <~ uj4t) <~ [~,(j),j = 1,2 . . . . . m;allt >~ t o
(2.1b)
ax(j) <~ xi(t) <<.[3x(i),i = 1,2 . . . . . n;allt >~ to,
(2.1c)
where f is a measured, b o u n d e d disturbance, x is an n-vector, u an m-vector, and ~ , F matrices of appropriate dimensions, c~and/~ define the controland state constraint sets. It is assumed that the whole state is measured (or estimated) at each sampling instant. More general linear time invariant constraints can also be introduced. The objective is to design a regulator that takes the state (or output) from its initial value to a prescribed stationary point, which will define the origin in the state space. It is assumed that the state and control constraints are admissible, i.e. the initial state is such that it is possible to stabilize (2.1) to the origin without violating the constraints. In G u t m a n (1984) necessary and sufficient conditions for admissibility are given when f = 0, and a variable structure linear state feedback controller is suggested. A straightforward way to solve the above control problem is to design a linear controller that never violates the constraints. Linear quadratic design may be used, or model following, i.e. the set point is changed in a modest way (e.g. a ramp reference input). Simulation must be used to check that these methods work. The c o m m o n method in practice is to design a linear controller disregarding the constraints and saturating the control input. The region of stability may be very small. For open loop stable plants Astr6m (1971) proposes a stabilizing relay, and Shapiro (1972) determines the matrix L so that u(t) = sat(Lx(t)) is stabilizing. Under certain conditions there exists a stabilizing u(t)= sat(Lx(t)) also for open loop unstable plants (Gutman, 1982a). Glattfelder (1974) and others have proposed the method of m a x - m i n selectors, whereby a set of controllers is used. W h e n a state constraint is violated, the control is switched to a controller
whose aim is to remove the state from the boundary. Although there is no proof of stability (the states may oscillate from one b o u n d a r y to another) the m e t h o d has become quite popular in e.g. the power industry. An elegant way to design a controller for (2.1) is to solve an optimal control problem with an optimization criterion reflecting the desire to stabilize the plant. In general, however, one gets an open loop control solution from one initial state, not the desired feedback control. To construct a feedback control, one may solve the optimal control problem for a great n u m b e r of initial conditions and then set up a look-up table. Alternatively, solve the optimal control problem on line at every sampling instant, and deliver the c o m m a n d signal pertaining to the current sampling interval.* This is the Open L o o p Optimal Feedback scheme (OLOF). Since the plant and constraints in (2.1) are linear, it is possible to use L P to solve the optimal control problem. However L P does not allow a free endtime. A free end-time problem can be solved by solving a sequence of fixed end-time problems, iterating over the end-time. Here the end-time is called "the time horizon r" and the end-time that solves the free end-time problem is called "the optimal time horizon Z'opt". Since it is assumed above that there exists a stabilizing control sequence, there exists at least one r for which the L P problem has a solution. So, at each sampling instant two or more optimization problems are solved for consecutive values of the time horizon z to find the predicted u(0), u(l) . . . . . u(z - 1), and x(1), x(2) . . . . . x(r) where r = 0 is present time. The optimization criterion is min ~. ( r l m a x [ - x i ( t )
+ 7.~(i),0]
i , t <<. z
+ wimax[xi(t ) - fix(i),0] + ~ cilxi('c)t,
(2.2)
i
where the vectors r and w are b o u n d a r y region penalty weights, ~'x and fix define the b o u n d a r y regions, and c is the terminal penalty weight. The optimal time horizon ropt is assumed to be found when Ix(r)la ~ ~', diixi(r)l ~< *:.
(2.3)
i=1
When solving the free end-time problem, it is essential that the iterations continue until topt is reached. A shorter final r might be harmful because stabilization to the origin might never be achieved. A longer final r can also be harmful in the O L O F setting, since the (non-unique) control obtained for the current sampling interval might counteract the control objective. However, during real time *The estimated computation time can be introduced as a constant delay in (2.1a).
A linear programming regulator operation it might happen occasionally that one sampling interval is not sufficient to solve the number of LP problems needed to find Topt. In such cases, the control suggested by the solution for the latest computed z is used. This is called "the suboptimal control". A standard linear programming package, PRIMAL, supplied by the Institute of Applied Mathematics, Stockholm (Holmberg, 1981) was used. The base inverse matrix is stored in product form, and reinversion is done when the multiplication becomes too time consuming, see Orchard-Hays (1968). The L P - O L O F controller was implemented on a VAX 11/780 computer, and the simulations done with the simulation language S I M N O N (Elmqvist, 1975). More details about the control algorithm are found in Gutman (1982a, b). 3. C O N T R O L
OF A LABORATORY
PROCESS
The LP regulator was used to control a laboratory process consisting of two water tanks with free outlets, as drawn in Fig. 1. The water heights were measured. The control objective was to get h2 to a reference value, h2ref , a s fast as possible, without causing overflow in either tank. The input voltage to the pump was limited. The non-linear pump characteristic was found to be so harmful that a local, analogue PID-feedback around the pump was introduced, whose dynamics were sufficiently fast to be neglected. Figure 2 depicts the block diagram of the process. The double tank system is non-linear. The time constant (from u to h2), when filling is about 125 s;
but only about 60 s when emptying. The sampling interval was set to 12s, because it is a suitable fraction of the shorter time constant, and because 12 s suffice as computation time for most reference changes. Dimensionless units were introduced: O~u~l,
~ Tank I
l
tocho signal I
= -A input voltagg
~
I
Pump [
]
Tonk
water FIG. I. The double tank process: physical set-up,
I
b2ref
FIG. 2. The double tank process; block diagram.
O~hl~l
,
O~h2~l,
(3.1)
u = 0 means no water flow, u = 1 means maximum flow. For ht and h 2 , 0 means that a tank is empty, 1 that it is full. The following continuous and sampled models were obtained from some step input experiments.
h(t) = \(--001970.0178 --0.01290) h(t) + _t/0"0263~_ LO )u(t) (3.2) "~ h(t + 12) = \0.175(0'7900.8570) h(t) + /'0.282 L0.0296/u(t) (3.3) where h(t) = [ht(t)h2(t)] v. This model was used in the LP regulator, with ULp ~ U - - Uref,
X 1 =- h l r e f ,
(3.4)
X 2 = h 2 -- h2ref
with hlref and u,,f computed from the stationary solution of (3.2), given h2ref. The input and state constraints are given in (3.1). In (2.2) and (2.3) e=0.05,
ftow
535
7x = [ - - h l r e f
r=w=[11], 0 ] T,
c=d=[l
6 x = [-hlref
Jr-
4],
0 . 9 0 ] x, (3.5)
were chosen. Since h2ref = 0 is an acceptable reference, which means ht = 0 in stationarity, there was no boundary penalty for low tank 1 levels. Therefore, y:,(1)= -hlref. The values of 7x(2) and 6x(2) mean that Ih2 - h2refl is included in the loss for all t, not only at the time horizon z. This is not the intended way to define a boundary region but it was decided on after a few preliminary experiments. The behaviour for small deviations from h2ref w a s especially improved. Moreover, when the sampling interval is not enough to find the optimal time horizon, it is good to penalize the deviation from h2ref along the trajectory to ensure that the control leads the system in the right direction. Since the second tank is essentially a first order system, this type of penalty function will not affect the time optimality of the LP solution. With 7x(2) 4:0 and 6,(2) :~ 0, the control tended to become a dead-beat control for small deviations from h2ref. This caused oscillations since the model
536
P.-O. GUTMAN
0
0
t000
0
1000
~___-f
.~__f
I reference
I/
i\ lo'oo
0
loICOmputing time [sees} in ir eoch sompling interval
O,L 0
,
10~)0
FIG. 3. The double tank process: u(t), hl(t ), hzref(t) ( = reference) and h2(t) during a real life experiment and the computation time in the simulated experiment.
was not perfect. It was also found that a minimum time horizon of two sampling intervals was beneficial for the local control, in addition to the above boundary region. The tanks were initially empty during the real life experiment, h2 was to follow h2ref , see Fig. 3. No outer disturbances affected the system. The computation times are perceptible in the figure. At reference value changes they were up to 4 s, except at the last reference value change, when the computations were interrupted after 10 s, and the suboptimal solution, u = 0, corresponding to • = 23, was used. The interruption was made because it was estimated that the next iteration would make the total computation time longer than 12 s. In "steady state", the computation times were less than 1.5 s. At reference value changes, the time horizon was restarted at 1. For reference value changes upwards, r roughly reflected the settling time, while it was too long for downward changes. This explains the overshoot at downward reference changes. The computation times at the reference changes were: for ~ = 4 , 0.23s; for • = 9 , 1.1s; for ~ = 11, 1.8 s; for ~ = 15, 3.4 s; and for r -- 23 (suboptimal), 10 s. These data suggest that the computation time is approximately proportional to r2 ln~. Notice that the state constraints are not violated, and that the local control is smooth. As a comparison, this LP regulator was simulated against the linear model (3.2). The computation times were included in the
simulation, but not in the sampled model. See Fig. 4. Notice the computation times in Fig. 3, which are roughly equal to the real life computation times. There are two main deviations between real life and simulation: in real life there is a slight stationary error in h2; the real life downwards dynamics are faster than the simulated. These discrepancies are due to the use of the simple linear model (3.3). The conclusions from the laboratory tank experiment, and supporting simulations, were, briefly: the L P - O L O F regulator can be used even if the computation time occasionally is longer than the sampling interval; the regulator functions well even if an approximate linear model of the non-linear
simulofed
I
/
/
h2(f)
//
0
0
1000
FIG. 4. The double tank process: h2(t) in the real life and simulated experiments.
A linear programming regulator plant is incorporated into it; the performance is sensitive to the tuning of the penalty function weights; the state constraints were rarely violated and the control constraints were never violated. The conclusions encouraged the author to find a real life process with suitable dynamics and essential control and state constraints on which the L P - O L O F regulator could be tested. The level control problem of the reservoir of the Bergeforsen hydroelectric station was thought to be suitable. The laboratory tank process could be seen as a very simple model of the reservoir. Therefore the approach to solve the Bergeforsen problem was similar: identification and modelling of a linear model, simulation, tuning of the penalty function weights and experimentation. 4. T H E B E R G E F O R S E N LEVEL C O N T R O L EXPERIMENT
The Bergeforsen reservoir, consisting of the lowest 55 km of Indals/ilven near Sundsvall, Sweden has at its upper end the Jiirkvissle power station and at its lower end the Bergeforsen power station. The flow at Bergeforsen ranges from ~ 5 0 m 3 s-1 to ~ 840m 3 s - i s and almost all the water enters from J/irkvissle. The Water Court decision, governing the operation of the power plants, states that the water level at K/ivsta, 14.3 km upstream from Bergeforsen "must coincide as close as possible but not exceed 23 m above sea level". See Fig. 5. At present the J/irkvissle power station is run according to the grid power demand while the four turbines at Bergeforsen are controlled manually to keep the K/ivsta level correct. Typically the J/irkvissle flow will be at a high level in the day and at a low level at night. The operators at Bergeforsen get oral information over the phone in advance of major changes in the J/irkvissle flow. The three upper curves of Fig. 7 represent normal operation over 75 h. A block diagram model of the automatic system is suggested in Fig. 6, with the J~irkvissle flow u 2 as a measured disturbance (feedforward), and the Kfivsta level y as the output feedback signal. The Bergeforsen flow Ul is the control signal. The control objective is to keep Jiirkvissle
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Kavsto level
Bergeforsen r---~----~ flow u~ ~ I Controller ~ Delay t - - ' - - ' ~ Dynamics I ~.___~r___a (m~/s) i , i i
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
y(mS/s)
.
.
.
.
.
.
.
.
J
FIG. 6. Block diagram of the Bergeforsen reservoir level control system.
y ~< 23 m, but as close as possible to 23 m in order to concur with the Water Court decision, and to utilize the reservoir capacity maximally. A transport delay is assumed from the inputs ul, u 2 followed by linear dynamics. Such a model is approximate since the reservoir dynamics are best described by non-linear partial differential equations. Previous data and physical insight suggest that the delay from ul should be 20min, the delay from u2 c a 1 h, and that the dynamics should include an integrator and some wave phenomenon. An identification experiment was performed, where especially u~ was varied randomly around the operating level, see Fig. 7. The sampling interval was set to 5 rain and necessary anti-aliasing prefiltering was done with analogue and digital filters. Correlation analysis, least squares and maximum likelihood methods of the interactive program Idpac (Wieslander, 1980) were used to identify the delays and dynamics. The identified model output is seen in Fig. 7. This model was modified when the L P - O L O F regulator was simulated against a computerized PDE model ot the reservoir, Routes, developed at the Swedish State Power Board. The final model was in state space form: x ( t + 1) = ~ x ( t ) + FUl(t ) + F(t)
(4.1a)
with x(t) = (y(t)ul(t ul(t -
- . 1)ul(t - 2)
3)ut(t - 4)ul(t - 5)
(4.1b)
u l ( t - 6)Aul(t - 2) Aul(t - 2)) r, zXul(t) = u l ( t + 1) -- u~(t).
,e..,[m]
< <
~ m31s Jarkvissle flow
.
Kavsta
[m3,s]
,,o.
537
5 5 Km
FIG. 5. The Bergeforscn reservoir.
Bercjeforsen
) )
(4.1c)
538
P.-O. GUTMAN
o
zoo
600
Somple 800
.oo
,oo
,oo
,oo
600
800
400
no,
,ooo
oJ[m]
ooT-- . . . . . . . . . . . . . . .
500 ; f ~ _ _ _ _ . v S ~
, Som~ no. [
o
,ooo
LK[ml
zz.91.-
_ '...
o
. .
Sample no.
~
20°
,o®
YM 0
.
0
.
200
.
.
400
600
l
800
I000
FIG. 7. Identification experiment, measurements for Bergeforsen flow (QB), J/irkvissle flow (QJ), KS.vsta level (LKI and idcntified model output for K/ivsta level with mean and linear trend subtracted (f'). 4.7451'10 -5
0
0
0-
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 0
0
1
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0 0
0 0
0 0
0 0
0 0
0 I
0 0 1
1
0
0
--5.7106" 10 -5
0
0
0
0
0
1
0
0
0
1
0
0
0
0
0
0
0
-1
0 0
-1
F(t)
= ((1.02-lO-Suz(t - 9)
1
-[-.f(uz,t
- 0 . 0 5 4 5 " 10 -5
i n p u t experiment) b/2o = 3 7 0 m 3 s - 1.
1)) 0 0 0 0 0 0 0 0)I",
+
(4.1e) f(u
2, t) =
lO-5(6.49(u2(t -
tl2(t -
15)) - 5.71(Uz(t - 15) - U z ( t -
- 0.78(uz(t - 16) --
3.8(U2(t
--
--
b~2(t
--
16))
(4.1f)
With the integrator emphasized, the model (4.1) is written y ( t + 1) = y ( t ) - 1.02" 10-5(ul(t -- 3)
17))
20) - u 2 ( t - 21))
+ 1.1(u2(t - 21) - u z ( t - 22)) + 2.7(Uz(t - 22) - u 2 ( t - 23))),
where u l , u2, y stand for Bergeforsen flow c o m m a n d input, JS.rkvissle flow a n d K&ivsta level a r o u n d their respective, n o m i n a l values. Ulo = 3 9 0 m 3 s - l Yo = 22.95 m (changed on line to 22.94 in step
(4.2)
The difference between Ulo a n d Uzo represents the average effect of side-flows into the reservoir. Aul was i n t r o d u c e d as two extra states in order to allow piecewise linear penalties on flow changes (see Fig. 8d). F ( t ) represents the feedforward term from u2, and f ( u 2 , t + 1) the wave generated by u2. Note that the s t a t i o n a r y value of f is 0.
14)
-
(4.1 d)
-
u2(t
-
9))
+ (terms representing waves) (see (4.1e) a n d note that
A linear p r o g r a m m i n g r e g u l a t o r
539
5
~b(l,i) = - 1 . 0 2 . 1 0 - 5 in (4.1d)).
The c o n s t r a i n t s a n d penalties are illustrated in Figs 8(a)-(d). Notice that high water levels are more penalized t h a n low (cf. W a t e r C o u r t decision), that the flow b o u n d s were c h a n g e d on line d e p e n d i n g on the n u m b e r of g e n e r a t o r s and that the introd u c t i o n of two Au t states m a d e it possible to shape the flow change p e n a l t y function, so that small changes are cheap, large changes costly a n d very large changes f o r b i d d e n (since they have d i s a s t r o u s effects on the delta d o w n s t r e a m from Bergeforsen). The m i n i m u m time h o r i z o n was set to 10, so that the current u 2 m e a s u r e m e n t w o u l d always enter (4.1), the m a x i m u m time horizon was set to 23, and the m a x i m u m C P U time before a s u b o p t i m a l c o n t r o l was used, 200s, The final c o n t r o l experiment e n c o m p a s s i n g 190 samples is shown in Fig. 9. T h e n u m b e r of g e n e r a t o r s was decided a c c o r d i n g to the p l a n n e d J/irkvissle flow, and the u~ b o u n d s were set accordingly. At 03.03 hours the Kfivstalevel reference Yo was c h a n g e d to 22.94 m. It is clear that y was c o n t r o l l e d satisfactorily, in general on the low side of Y0, with an acceptable oscillative a m p l i t u d e (cf. Fig. 7 when the control was manual).
i=3
T o tune the p a r a m e t e r s of the L P - O L O F controller, it was f o u n d to be beneficial to simulate it against m o d e l (4.1), i.e. the m o d e l internal to the controller. This is called i n t e r n a l s i m u l a t i o n . M o r e over this m o d e l can be used to p r e d i c t future o u t p u t s d u r i n g o p e r a t i o n . W h e n the p r e d i c t i o n time o r the time h o r i z o n z of the c o n t r o l a l g o r i t h m is longer t h a n 10 s a m p l i n g periods, then future, as yet u n m e a s u r e d u 2 values enter (4.1). u2 was then predicted in the simplest possible way: U2('C + g) = Uz(t), Z > 0, where u2(t ) = latest m e a s u r e m e n t of u 2. In a refined controller, future k n o w n changes of u 2 can be a c c o u n t e d for. T h e final p a r a m e t e r set was as shown in Table 1, a n d % = Ctx(i ),
i = 1,2 . . . . . 7
ft, = fix(i),
i = 1,2 . . . . . 7.
(4.3)
O n the first occasion, when 3 turbines were used, the u p p e r u l - b o u n d was 230, the second time it was 240. It can be shown that for limited u2 a n d Au 2, (4.2) a n d T a b l e 1 are admissible c o n s t r a i n t sets. In (2.2), (2.3), ci
=
10;
d i = ci;
rt
=
200;
to i = 1000;
c i = 0,
i=
i = 2 . . . . . 9;
1. . . . . 9;
Ti = 0 ,
i = 2 . . . . . 7;
O,)i ~ ri~
i = 2 . . . . . 9;
r 8 = 0.1,
r9
=
0.01;
(4.4)
= 0.05. TABLE 1.
S t a t e no. (5.1)
~tx(i) (2.1c)
)'x(i) (4.2)
6~(i)
flxll)
(4.2)
(2.1c)
Level Flow
1 2..... 7
zXu* Au,*
8 9
-0.15 - 320 - 290 -100
-0.004 - 320 - 290 -15 0
0.004 20 230/240 15 0
0.15 20 230/240 100 100
Eqn. no:
* nominal
value
-100 =
Operating condition: no. of turbines 2 3
0
Intermediate
Terminal level penalty
level
penalty
10o
I
' I
-0.15
stopping criterion
[m] 0.15
toleronee
FI(;. 8(a), P e n a l t y f u n c t i o n c o m p o n e n t s : t e r m i n a l level penalty.
I
-0.15
I
0
QI5
FIG. 8(b). Penalty function components: intermediate level penalty.
;-
540
P.-O. GUTMAN One exception is the overshoot at 08.45 of 3.4cm which did not fall within the unofficial tolerance at 3 cm. It is believed that the reason is the insufficiency of the model; in particular, the time delay from u~ should have been 35 min instead of 20 min. The mismodelling is probably due to weak identifiability at the delay. In the identification experiment (Fig. 7) u~ changed only a couple of times. In the judgement of the operators who run the plant manually, the L P - O L O F control was satisfactory, corresponding to an average operator. The advantages were seen as (1) the possibility to change Bergeforsen flow bounds on line; (2) the flow bounds were not violated; (3) the prediction of future K/ivsta levels and control inputs, the latter making it possible to govern the flow (suboptimally) for a few sampling intervals even when the computer went down (this happened at one occasion); (4) the feedforward from the J/irkvissle flow; (5) the nonsymmetrical penalty on y and the possibility to restrict flow changes. All in all, with somewhat more exact modelling the L P - O L O F regulator seems viable for the reservoir level control problem.
Flow penalty '(2 turbine case)
I -300
I -200
I -I00
I I00
)
FI(,. 8(c). Penalty function components: flow penalty. Flow c h a n g e 'penalty
l
,o
5. C O N C L U S I O N S
Of course, few general conclusions about the LPO L a F regulator can be drawn from two experiments. The most obvious advantage over most alternative methods is that the state and control constraints are taken into account when the control is computed, and that a criterion reflecting the design goal is optimized. Some other general points arc discussed in G u t m a n 1982, 1983). Here some specilic remarks are made concerning the on-line
3/S/5rain] I
-Io0
o
)
I0o
I"1(;. 8(d). Penally, function components: flow change penalty.
,[m3/s] f l o w
r
600 500. . . . "~"
QB-bounds
400-
J
. . . . . . . . .
....... ........
500200I00time[hours]
t 02.15
I
2015 ,[m]
2215
0015
I
0415
I 0915
I
0615
I
level
2297
:l :' "'-..
22 95
:
,'.. ", ......: ......L K
,-.Y"., \,.
\
22 95-
,/.Y-
_
:""... / ...... • :
: ."
I
:, L K r e f ., -"...................
..
. :,.:
""":; 2015
)
1015
2215
O015
L,
'." A .........-
............
time[hours] 02.15
0415
06.15
0815
10.15
I't~;. t.L l)iagran~ of Imal cxperimcm: Bergcforsen I1~',', ~QB), Jiirk,,isslc flow (QJ). Kfivsta level H_KI and Kii,,sla rclk:rcnce level (kKrcf).
A linear p r o g r a m m i n g r e g u l a t o r e x p e r i m e n t s . T h e r e g u l a t o r p r o v e d itself r o b u s t f r o m t h r e e p o i n t s o f view. (1) I n spite o f the a s s u m e d h i g h l y a p p r o x i m a t e d l i n e a r m o d e l s of the n o n l i n e a r plants, it p e r f o r m e d well. (2) It w i t h s t o o d well o c c a s i o n s w h e n the c o m p u t a t i o n t i m e s were l o n g in c o m p a r i s o n to the s a m p l i n g interval, by u s i n g s u b o p t i m a l c o n t r o l s ; a n d a n o c c a s i o n of computer shutdown when previously computed p r e d i c t e d c o n t r o l s w e r e used. (3) It is sensitive to t h e t u n i n g of the p e n a l t y f u n c t i o n weights. T h e B e r g e f o r s e n e x p e r i m e n t s h o w e d t h a t o n e of the m a i n a d v a n t a g e s of the L P - O L O F r e g u l a t o r is the p o s s i b i l i t y o f h a v i n g h a r d b o u n d s o n the c o n t r o l variables and their derivatives, which are taken into a c c o u n t w h e n the c o n t r o l a c t i o n is c o m p u t e d (cf. w i t h s a t u r a t o r s a c t i n g o n the a l r e a d y c o m p u t e d control). O t h e r advantages include dynamic c h a n g e s of b o u n d s , s k e w e d p e n a l t y functions, pred i c t i o n s o f f u t u r e p r o c e s s values a n d the possibility of t u n i n g the L P - O L O F by s i m u l a t i n g it a g a i n s t itself. Acknowledgements--The Bergeforsen experiment was performed
under the aegis of Kockumation AB, MalmS, (thanks to J. Tuszynski) and the Swedish State Power Board (thanks to A. Bruse and G. Fabricius). The report was completed during the author's stay at the Technion-Israel Institute of Technology, which was made possible by generous grants from the C. F. LundstrSm Foundation, the Swedish Board for Technical Development, the Helge Ax:son Johnson Foundation and the Per Westling Foundation of Lund University. The author thanks the reviewers, in particular Reviewer # 2, who suggested several relevant references on water systems control.
7. REFERENCES AstrSm, K. J. (1971). Olinj~ira system (in Swedish). Lund Inst. Techn., TLTH/VBV, Lund, Sweden. Bashein, G. (1971). A simplex algorithm for on-line computation of time optimal controls, lEE Trans. Aut. Control, AC-16, 479-482. Cadzow, J. A. (1974). Minimum amplitude control of linear discrete systems. Int. J. Control, 19, 765. Canon, M. D., C. D. Cullum and E. Polak (1970). Theory of Optimal Control and Mathematical Programming. McGrawHill, New York. Chang, T. S. and D. E. Seborg (1983). A linear programming approach for multivariable feedback control with inequality constraints. Int. J. Control, 37, 583-597. Dreyfus, S. (1964). Some types of optimal control of stochastic systems. S l A M J. Control, 2, 120. Elmqvist, H. (1975). Simnon An interactive simulation program for nonlinear systems. Rep. 7502, Dept. Aut. Control, Lund Univ., Sweden. Glattfelder, A. H. (1974). Regelungsysteme mit Bergrenzungen. Oldenburg, Mi.inchen. Gutman, P.-O. (1982a). Controllers for bilinear and constrained linear systems. Ph.D. dissertation, Dept. Aut. Control, Lund Univ., Sweden, CODEN:LUTFD2/TFRT- 1022)/1-133/ (1982).
AUTO 22:5-C
541
Gutman, P.-O. (1982b). On-line use of a linear programming controller. 3rd IFAC/IFIP Symp. Software Comp. Control, Madrid, Spain, October 5-8. Gutman, P.-O. (1983a). On the linear programming OLOFregulator. Technion-Israel Inst. Tech., Haifa, Israel, EE Pub. No. 469. Gutmam P.-O. (1983b). Test of the OLOF-regulator for water reservoir level control. Technion-Israel Inst. Tech., Haifa, Israel, EE Pub. No. 468. Gutman. P.-O. (1984). Admissible sets and feedback control for discrete time linear dynamical systems bounded controls and states. Proc. 23rd IEEE Conf. Decis. Control, December 1214, Las Vegas, Vol. 3, pp. 1727-1731; IEEE l~'ans. Aut. Control, AC-31, pp. 373-376 (1986). Holmberg, G. (1981). PRIM2-ett subrutinpaket som 15ser LPproblem (in Swedish). UTM-Prog. nr. 50. Uppdragsgruppen f. tillfimpad matematik. Royal Inst. of Techn., Stockholm, Sweden. Kolev, L. V. (1978). Minimum-fuel minimum-time control of linear discrete systems. Int. J. Control, 27, 21-29. Knudsen, J. K. H. (1975). Time-optimal computer control of a pilot plant evaporator. Proc. IFAC 6th World Congress, Part 2, pp. 45.3/1-6. Larson, R. E. and W. G. Keckler (1969). Application of dynamic programming to the control of water resource systems. Automatica, 5, 15-26. Loucks, D. P. (1968). Computer models for reservoir regulation. J. Sanit. Eng. Div., Amer. Soc. Civil Eng., 94(SA-4), 657-669. M aas, A. et al. (1962). Design of water resource systems. Harvard University Press, Cambridge, Massachusetts. Manne, A. S. (1960). Linear programming and sequential decisions. Manag. Sci., 6, 259-267. Nayak, S. C. and S. R. Arora (1971). Optimal capacities for a multi-reservoir system using the linear decision rule. Water Resource Res., 7, 485 498. Orchard-Hays, W. (1968). Advanced Linear Programming Computing Techniques. McGraw-Hilk New YorL Propoi, A. I. (1963). Use of L P methods for synthesizing sampleddata automatic systems. Aut. Remote Control, 24, 837. Rasmy, M. E. and M. H. K. Hamza. Minimum-effort timeoptimal control of linear discrete systems. Int. J. Control, 21, 293-304. Revelle, C., D. P. Loucks and W. R. Lynn (1968). Linear programming applied to water quality management. Water Resource Res., 4, 1-10. Revelle, C., E. Joeres and W. Kirby (1969). The linear decision rule in reservoir management and design I. Development of the stochastic model. Water Resource Res., 5, 767-777. Sachdeva, S. A. (1982). Bibliography on optimal reservoir drawdown for the hydroelectric-thermal power system operation. IEEE Trans. Power Apparat. Syst., PAS-101, 14871495. Shapiro, E. Y. (1972). Suboptimal design of minimal-time feedback regulator for linear time-invariant plants. Sc.D. Dissertation, Dept. Elec. Eng. Comp. Sci., Columbia Univ., New York. Vlieger, J. H. de, H. B. Verbruggen, P. M. Bruijn (1981). A time-optimal control algorithm for on-line computer control allowing bounding control and state variables. Control Comp., 9, 1-6. Vlieger, J. H. de, H. B. Verbruggen, P. M. Bruijn (1982). A timeOptimal control algorithm for digital computer control. Automatica, 18, 239-244. Wieslander, J. (1980). Idpac commands--user's guide. Dept. Aut. Control, Lund, Sweden, CODEN: LUTFDZ/(TFRT3157)/1-108/(1980). Zadeh, L. A. and B. H. Whalen (1962). On optimal control and linear programming. IRE Trans. Aut. Control, 7, 45.