- Email: [email protected]
A Moving Horizon Approach to a Chance Constrained Nonlinear Dynamic Optimization Problem*
A Moving Horizon Approach to a Chance Constrained Nonlinear Dynamic Optimization Problem ⋆ Abebe Geletu*, Michael Kl¨ oppel, ∗ Pu Li*, Armin Hoffmann ∗∗ ∗ Ilmenau University of Technology, Institute of Automation & Systems Engineering, Postfach 10 05 65, 98684 Ilmenau, Germany. (e-mail: [email protected]) (e-mail: [email protected] ) (e-mail: [email protected]) ∗∗ Ilmenau University of Technology, Institute of Mathematics, Dept. of OR and Stochastic, Postfach 10 05 65, 98684 Ilmenau, Germany. (e-mail: [email protected]).
Abstract: In this paper, a chance constrained nonlinear dynamic optimization problem is considered, which will be investigated by using a moving horizon scheme. In each horizon, the chance constraints will be written (transformed) in terms of those (input) random variables with known probability distributions by using monotonicity relations. Some definitions and properties related to the required monotonicity properties are introduced. For the application problem considered these monotonicity properties hold automatically true. The chance constraints and their gradients are evaluated by computing multivariate normal integrals using direct numerical integration. Numerical experimentation results will also be reported. Keywords: Chance Constrained Optimization, Dynamic Optimization, SQP, MPC, Montonic Relation, Normal Distribution, Numerical Integration. 1. INTRODUCTION In this paper the following nonlinear optimization problem is considered min
f (u)
(1)
s.t. g(u, y, ξ) = 0 P r{G(u, y, ξ) ≤ 0} ≥ α, α ∈ [0, 1],
(2)
straints. Currently these models are being widely studied and applied in chemical process engineering (Li et al. (2002), Li (2007), Li et al. (2008), Bartl (2006), Flemming et al. (2007), Henrion et al. (2001), Wendt et al. (2002)); optimal power flow (Shiina (1999)), finance and economics (Chanres and Cooper (1959)), water resource management (Prekopa (1995)), reliability based design optimization (Chan et al. (2006)), etc. Special cases (q2 = m) of (4) and (5) are either
(3)
where u, ξ, y are the vectors of the decision, uncertain input and uncertain output variables, respectively. The probability P r(·) is a measure of the confidence level of satisfying the inequality constraints; where α ∈ [0, 1] is a pre-given confidence (reliability) level. Stochastic optimization problems of this form are commonly known as chance (or probabilistic) constrained optimization problems. Assuming that f : Rn → R, g : Rn × Rm × Rl → Rq1 , and G : Rn × Rm × Rl → Rq2 , the chance constraint can have one of the following two forms: either P r{Gi (u, y, ξ) ≤ 0} ≥ αi , i = 1, . . . , q2 , αi ∈ [0, 1] (4) or P r{Gi (u, y, ξ) ≤ 0, i = 1, . . . , q2 } ≥ α, α ∈ [0, 1]. (5) The former are called separate or single chance constraints; while the latter are known as joint chance con⋆ We are very grateful to the DFG - Deutsche Forschungsgemeinschaft for supporting this research with Project No. LI806/8-2.
P r{yimin ≤ yi (u, ξ) ≤ yimax } ≥ αi , i = 1, . . . , m,
(6)
or P r{yimin ≤ yi (u, ξ) ≤ yimax , i = 1, . . . , m} ≥ α. (7) These constraints define the level of confidence that bounding constraints should be satisfied for a reliable as well as optimal operation of a system. The numerical computation of such chance constrained optimization problems is quite difficult and expensive. The usual approach is by relaxation of the problem into a deterministic non-linear optimization problem. 2. MONOTONICITY APPROACH TO THE SOLUTION OF CHANCE CONSTRAINED OPTIMIZATION PROBLEMS Of interest is the chance constrained nonlinear optimization problem (1) - (3) with either of the special chance constraints (6) or (7). It is assumed that the uncertain
variables ξ posses some joint probability density ρ(ξ). Such a density could be determined based on historical data of the application problem being considered. However, due to the (possible) nonlinear relationships between the variable y and ξ, it is usually hard to find out the probability density of the output variables y. Thus, the idea proposed by Wendt et al. (2002) is to transform the chance constraints in the space of y into chance constraints in the space of the variable ξ whose joint distribution is already known. The probability-measure preserving transformation can be guaranteed through a monotonic relation among the uncertain input variables ξ and output variables y. In case of separate chance constraints, it is assumed that for each of the restricted variables yi , there is a function F relating yi to one of the inputs ξj(i) such that ξj(i) 7→ fi (ξj(i) ) := F (ξ1 , . . . , ξj(i) , . . . , ξl , u) = yi
That is, there is a monotonic relation between yi and ξj(i) ; i.e., as ξj(i) increases, increases yi , too; ( it can also be the case that as ξj(i) increases, decreases yi ). There are several applied optimization models in which such monotonic relations readily available. (See the model in Section 4). In fact, the approach is to use the equality constraints g(u, y, ξ) = 0 in verifying such monotonicity relations between the variables yi and ξi ; furthermore, for a given u, to solve the systems of equations g(u, y max , ξ) = 0,
in order to determine ξ min (u) and ξ max (u), respectively; so that the (say, separate) chance constraints can be equivalently written as P r{yimin ≤ yi (u, ξ) ≤ yimax } min max = P r{ξj(i) (u) ≤ ξj(i) ≤ ξj(i) (u)} ≥ αi ,
for each i = 1, . . . , m. By doing so the equality constraints g(u, y, ξ) = 0 will be effectively eliminated from the optimization problem and the separate chance constrained problem will be transformed into min s.t. Z
f (u) +∞
...
Z
(8)
max ξj(i) (u)
...
min (u) ξj(i)
−∞
Z
+∞
ρ(ξ1 , . . . , ξl )dξl . . . dξ1 (9)
By defining hi (u) :=
+∞
... −∞
∇hi (uk ) =
Z
+∞
−∞ max k ∂ξj(i) (u )
...
Z
+∞
−∞
−
h
max k ρ(ξ1 , . . . , ξj(i) (u ), . . . ξl )×
min k ρ(ξ1 , . . . , ξj(i) (u ), . . . ξl )
min k ∂ξj(i) (u )
∂u
#
max ∂ξj(i) (uk ) ∂u
∂ξ min (uk )
and j(i) can be respecThe derivatives ∂u tively determined from the equation g(u, ξ, u) = 0 by the implicit function theorem. Consequently, the computation of h(uk ) and ∇hi (uk ) involves the evaluation of integrals. In fact, this approach requires efficient tools for numerical computation of multi-dimensional integrals. Properly speaking, the efficiency of the approach rests upon the precise computation of integrals. For this reason quadrature based methods, such as collocation on finite elements, are quite preferred instead of methods based on sampling. 3. DIRECTIONAL MONOTONICITY The above solution approach relies on a monotonic relationship between the input and output variables. Therefore, it is required to conduct a monotonicity analysis of a given system. However, such monotonic relations may not be available in all cases. Recently, monotonic properties of continuous dynamical systems are studied by Angeli and Sontag (2003). They give a general monotonicity property based on tangent cones; which they later on specialize by using partial derivatives. Here is given a similar, but weaker monotonicity property. Now, considering chance constraints (6) and for the sake of the practical problem to be considered in the next section, it is assumed that the constraint g(u, y, ξ) = 0 to have the form yj = gj (u, ξ1 , . . . , ξl ), j = 1, . . . , m; (i.e. m = q1 = q2 ).
−∞
≥ αi , i = 1, . . . , m, αi ∈ [0, 1].
Z
hi (u) ≥ 0, i = 1, . . . , m. This problem can be non-convex. Thus, an efficient numerical method is required for its solution. Currently, the method of choice for such nonlinear programming problem (NLP) is the interior-point method (see W¨achter and Biegler (2006))). However, to apply this NLP algorithm, at each step k, given uk , the constraint values hi (uk ) and gradients ∇hi (uk ) must be computed, where (see Xie et al. (2007), also Prekopa (1995))
dξl . . . dξj(i)−1 dξj(i)+1 . . . dξ1 .
ξj(i) = fi−1 (yi ).
and
f (u)
s.t.
∂u
and fi (·) is strictly monotonic; hence,
g(u, y min , ξ) = 0
min
Z
max ξj(i) (u)
...
min (u) ξj(i)
Z
+∞
ρ(ξ1 , . . . , ξl )dξl −∞
. . . dξ1 ) − αi , i = 1, . . . , m,
the problem (8)-(9) is written as
Definition 1. (Directional Monotonicity on a segment). Let g : Rn → R, y = g(ξ1 , ξ2 , . . . , ξn ) and a, b ∈ Rn , such that b > a . A function g is said to be monotonically increasing on the segment [a, b] in the direction of v = b−a iff the function ϕ(λ) = g(a + λv) is montonically increasing on [0, 1]. Theorem 2. Suppose g : Rn → R is a differentiable function and v ∈ Rn . Then g is monotonic increasing on the segment [a, b], a < b, if and only if ∇g(ξ)⊤ v ≥ 0, ∀ξ ∈ [a, b].
The proof follows trivially from the well known mean-value theorem (see for instance Spivak (1994)). The above directional monotony can be specialized to monotonicity in coordinate directions. In this case, let v be one of the standard unit vectors; i.e. v = ei (i ∈ {1, . . . , n}). Thus, ϕ(λ) = g(ξ + λei ) = F (ξ1 , . . . , ξi + λ, . . . , ξn ). As a result, g is monotonic on the segment [a, b] in the direction ei if and only if ϕ′ (λ) = λ∇g(ξ)⊤ ei ≥ 0, λ ∈ [0, 1]. Consequently, g is increasing on the segment [a, b] in the direction of the unit vector v = ei iff ∂g(ξ) ≥ 0, ∀ξ ∈ [a, b]. ∂ξi Accordingly, for y = g(ξ1 , ξ2 , . . . , ξn ), if it holds that ∂g(ξ) ∂ξi ≥ 0, then the variable y increases as the ith coordinate ξi increases. If y, ξ1 , . . . , ξn are random variables with some probability distributions, then it follows that P r {aj ≤ yj ≤ bj } = P r g −1 (ξ1 , . . . , ξi−1 , aj , ξi+1 , . . . , ξn ) ≤ ξi ≤ g −1 (ξ1 , . . . , ξi−1 , bj , ξi+1 , . . . , ξn )
The monotonicity properties being used in the next section can be interpreted as given above. Furthermore, with respect to the joint chance constraint (7), a monotonic relation should guarantee the transformation of a box constraint in the space of output variables y into a box constraint in the space of corresponding uncertain input variables ξ. In general, this is not a trivial issue. Considering the existing non-linear relation between y and ξ, work is still in progress to find out a viable monotonicity property.
Fig. 1. Buffer tank
a finite number of unit time steps t = 0, 1, . . . , N − 1, this problem is described as follows:
Objective f unction : φ(u0 , ..., uN −1 ) → min State equations : t = 1, . . . , N ; Vt = Vt−1 + qt−1 − ut−1 ; qt−1 Ct = Ct−1 + [COt−1 − Ct−1 ] ; Vt Initial conditions : V0 = v0 , C0 = c0 Control constraints : t = 1, . . . , N ; umin ≤ ut ≤ umax ;
4. AN OPTIMAL CONTROL PROBLEM WITH CHANCE CONSTRAINTS In a chemical process the feed to a distillation column can be an outflow of several upstream plants. For a smooth operation of the distillation column it is common to use a buffer tank to guarantee continuous supply to the column and avoid process disruption as shown in Fig. 1. Due to upstream processes there might be (uncertain) inflow disturbances to the buffer tank. The accumulation in the buffer tank should have certain minimum level (to guarantee continuity of the distillation process) and should not exceed the tank capacity. Moreover, the content of the buffer tank will be regulated by the outflow; i.e. the inflow to the distillation column. Let qt = q(t) represent the uncertain feed flow rate and COt = CO(t) inflow component concentration into the tank. The distillation column takes, at a time instant t, an amount ut = u(t) from the buffer tank for further processing. The volume Vt = V (t) and the concentration Ct = C(t) in the tank should satisfy bounding constraints up to a certain pre-given confidence levels α1 and α2 , respectively. The initial state of the tank, at t = 0, is given by V0 = V (0) and C0 = C(0). It is desired that the amount taken from the tank, from one time point to the next, must have a minimum fluctuation. Considering
Chance constraints : t = 1, . . . , N ; P {Vmin ≤ Vt ≤ Vmax } ≥ α1 ; P {Cmin ≤ Ct ≤ Cmax } ≥ α2 . Here, it is assumed that the uncertain inputs to have normal distributions qt ∼ N (µ(qt ), σt2 ) and COt ∼ N (µ(COt ), ̺2t ), where σ 2 and ̺2t are given variance of qt and COt , respectively, and µ(·) represents expected values. Linear problems of the above type have been considered, for instance, by Geibel and Wysotzki (2005), Henrion et al. (2001), Li (2007); Li et al. (2002); non-linear problems by Li et al. (2008) and Xie et al. (2007). In Li et al. (2002, 2008) a truncated normal distribution (see Horace (2005)) for the uncertain input variables has been used, while Henrion et al. (2001) studied the problem by using a scenario approach. In this paper, the uncertain variables are assumed to have standard normal distributions without truncation and integrals are evaluated by the direct numerical method. The approach used to solve the above problem is the moving horizon scheme. The control and prediction horizons are taken to be of equal length h. Thus, for each k = 0, . . . , N − 1, the following problem is to be solved:
min (∆u)⊤ ∆u s.t. i−1 X
56
(qk+j − uk+j ), i = 1, . . . , h
j=0
Ck+i
V − u k+j k+j Ck = V k+j+1 j=0 i−1 i−1 X Y Vk+s − uk+s qk+j COk+j , + V Vk+j+1 k+s+1 j=0 s=j+1
i−1 Y
54 C0(k) in g/l
Vk+i = Vk +
Different realizations of the concentration C0(k)
58
52 50 48 46 44
1
2
3
4
5
6
7
8
9 10 11 Time step k
12
13
14
15
16
17
18
19
i = 1, . . . , h;
umin ≤ uk+i ≤ umax , i = 1, . . . , h; P r{Vmin ≤ Vk+i ≤ Vmax } ≥ α1 , i = 1, . . . , h;
Fig. 2. Realization of the uncertain inflow concentration of COt
P r{Cmin ≤ Ck+i ≤ Cmax } ≥ α2 , i = 1, . . . , h;
=⇒ Vk+i ↑
and COk+i−1 ↑
=⇒ Ck+i ↑,
i = 1, . . . , h; meaning that Vt increases as qt−1 increases or Vt decreases if qt−1 decreases, etc. When transforming the chance constraint for Vk+i , the uncertain variables qk , . . . , qk+i−1 are taken into consideration; while w.r.t. Ck+i , the variables qk ,. . .,qk+i−1 , COk ,. . .,COk+i−1 are considered. Thus, in each horizon, it is assumed that each group of stochastic variables qk ,. . ., qk+i−1 and qk ,. . ., qk+i−1 , COk ,. . .,COk+i−1 have their respective multivariate normal distributions with corresponding covariance matrices. Accordingly, the length of the horizon determines the number and dimension of integrals to be computed. Furthermore, derivative values of constraints, required for the NLP algorithm, can be evaluated either iteratively using the method suggested by Prekopa for multivariate normal integrals (see Prekopa (1995) ) or by using quadrature methods (Bartl (2006)). In this paper a particular case of h = 3 for the purpose of numerical experimentation is considered. The following data are defined: • initial values V0 = 160, C0 = 50; • bounds Vmin = 130, Vmax = 170, Cmin = 49, Cmax = 51, umin = 0, umax = +∞; • reliability levels α1 = α2 = 0.8; • µ(qt ), µ(COt ), σt = V ar(qt ) = 0.7, ̺t = V ar(COt ) = 1, t = 1, . . . , N • N = 20. In each horizon the non-linear optimization problem is solved using IPOPT (W¨ achter and Biegler (2006)). Values of chance constraints and their gradients are computed using codes developed by Bartl (2006) based on orthogonal collocation. Integration limits are determined by solving the equality constraints using the Newton method. The following figures depict the results of numerical experiments for 21 realizations. Comparing Fig. 2 and Fig. 3, the realizations for COt and Ct justify the monotonic relations between them. It can
50.6 50.4 C(k) in g/l
qk+i−1 ↑
50.8
50.2 50 49.8 49.6 49.4 49.2 49
2
4
6
8
10 12 Time step k
14
16
18
20
Fig. 3. Resulting tank concentration for the 20 time intervals Different realizations of the feed q(k)
18 16 14 12 Feed q(k) in l
where ∆u = (uk+1 − uk , . . . , uk+h − uk+h−1 ), and define Qt−1 s=t (·) = 1. It is easy to observe that the following monotonic relationships hold true:
Resulting tank concentration C(k) for realizations of q(k) and C0(k) and calculated u(k)
51
10 8 6 4 2 0
1
2
3
4
5
6
7
8
9 10 11 Time step k
12
13
14
15
16
17
18
19
Fig. 4. Realization of the uncertain feed flow qt also be seen from Fig. 4 and Fig. 5 that the volume Vt of the buffer tank relates with the inflow qt monotonically. It is not surprising to observe from Fig. 5 and Fig. 6 that the control ut imitating the behavior of the tank volume Vt : when there is little in the tank, little will be allowed to flow out; and when there is more, then more will be allowed to flow out. The realized values of tank volume Vt (Fig. 5) and concentration Ct (Fig. 3) satisfy, in almost all the 21 realizations, their bounding constraints up to the required reliability levels α1 = 0.8 and α2 = 0.8. Good choice for α1 and α2 are values greater or equal to 0.5 (at least 50% reliability). Surely a process is expected to work with a higher relia-
REFERENCES
Resulting tank volume V(k) for realized q(k) and calculated u(k)
180
170
V(k) in l
160
150
140
130
120
2
4
6
8
10 12 Time step k
14
16
18
20
Fig. 5. Realized tank volume Resulting controls u(k) for the realizations of q(k) and C0(k)
18 16 14
u(k) in l
12 10 8 6 4 2
1
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16 17 18 19 20 Time step k
Fig. 6. Control realization for the 20 time intervals bility level; i.e. with values for α1 and α2 near 1. However, requiring higher reliability levels may cause infeasibility of the constraints. The analysis of the largest values of α1 and α2 that still guarantee the feasibility of the problem has not been done yet. This will be a future undertaking. A similar investigation for linear systems can be found in Li et al. (2002). The numerical experiments in this work were carried out using - Matlabr R2008a for Linux 64-bit, - IPOpt Version 3.4.2 with the linear solver ma27, - Processor: Athlon64 X2 3800+, 2GB RAM. Integrals are calculated with a maximum absolute error of 0.002 and the NLP algorithm is allowed to work with a termination tolerance of 0.0001. Individual realizations required an average total computing time of approximately 100 minutes. Among these, most of the time spent went for the calculation of integrals. CONCLUSIONS: Nonlinear dynamic optimization under chance constraints is considered. The solution of the problem depends on a monotonicity relation between uncertain input and output variables. For this purpose a monotonicity analysis is introduced to transform chance constraints from the space of uncertain output variables into the space of uncertain input variables. The moving horizon scheme is used for the optimal control of a tank process under the feed flow and concentration uncertainties.
Angeli, D. and Sontag, D. (2003). Montone control systems. IEEE Trans. on Auto. Ctrl., 48, 1684–1698. Bartl, M. (2006). Station¨are Prozessoptimierung unter Einbezug von unsicheren Randbedingungen am Beispiel einer Destillationsanlage. Diplomarbeit, Fakult¨ at f¨ ur Informatik und Automatisierung, Technische Unviersit¨ at Ilmenau. Chan, K.Y., Skerlos, S.J., and Papalambros, P.Y. (2006). Optimal design with non-normally distributed parameters, conditional probability, and joint constrained reliabilities. Proceedings of the 2006 ASME International Design Engineering Technical Conferences and Computers and Information In Engineering Conference, Philadelphia, PA, DETC2006/DAC-99417, 14 Pages. Chanres, A. and Cooper, W. (1959). Chance-constrained programming. Management Science, 6, 73–79. Flemming, T., Bartl, M., and Li, P. (2007). Set-point optimization for closed-loop control systems under uncertainty. Industrial & Engineering Chemistry Research, 46, 4930–4942. Geibel, P. and Wysotzki, F. (2005). Risk-sensetive reinforcement learning apllied to control under constraints. J. Artificial Inteligence Res., 24, 81–108. Henrion, R., Li, P., M¨oller, A., Seinbach, M.C., Wendt, M., and Wonzy, G. (2001). Stochastic optimization for operating chemical process under uncertainity. In M. Gr¨othschel, S. Krumke, and J. Rambau (eds.), Online Optimization of Large Scale Systems, 455 – 476. Springer Verlag. Horace, W.C. (2005). Some results on the multivariate truncated normal distribution. Journal of Multivariate Analysis, 94, 209–221. Li, P. (2007). Prozessoptimierung unter Unsicherheiten. Oldenbourg Verlag. Li, P., Arellano-Garcia, H., and Wozny, H. (2008). Chance constrained programming approach to process optimization under uncertainity. Computers and Chemical Engineering, 32, 24–45. Li, P., Wendt, M., Arellano-Garcia, H., and Wozny, G. (2002). Optimal operation of distillation process under uncertain inflows accumulated in a feed tank. AIChe Journal, 48, 1198–1211. Prekopa, A. (1995). Stochastic programming. Kluwer Academic Publishers. Shiina, T. (1999). Numerical solution technique for joint chance-constrained programming problem- application to electric power capacity expansion. J. of OR Socity of Japan, 42, 128 – 140. Spivak, M. (1994). Calculus, 3rd ed.,. Publish or Perish Inc. W¨achter, A. and Biegler, L.T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Prog., 106, 25 – 57. Wendt, M., Li, P., and Wozny, G. (2002). Nonlinear chance-constrained process optimization under uncertainity. Ind. Eng. Chem. Res., 41, 3621 – 3629. Xie, L., Li, P., and Wozny, G. (2007). Chance constrained nonlinear model predictive control. In R. Findeisen, F. Allg¨ower, and L.T. Biegler (eds.), Assessment and Future Directions of nonlinear model predictive control, 295 – 304. Springer Verlag, London.
Abstract: In this paper, a chance constrained nonlinear dynamic optimization problem is considered, which will be investigated by using a moving horizon scheme. In each horizon, the chance constraints will be written (transformed) in terms of those (input) random variables with known probability distributions by using monotonicity relations. Some definitions and properties related to the required monotonicity properties are introduced. For the application problem considered these monotonicity properties hold automatically true. The chance constraints and their gradients are evaluated by computing multivariate normal integrals using direct numerical integration. Numerical experimentation results will also be reported. Keywords: Chance Constrained Optimization, Dynamic Optimization, SQP, MPC, Montonic Relation, Normal Distribution, Numerical Integration. 1. INTRODUCTION In this paper the following nonlinear optimization problem is considered min
f (u)
(1)
s.t. g(u, y, ξ) = 0 P r{G(u, y, ξ) ≤ 0} ≥ α, α ∈ [0, 1],
(2)
straints. Currently these models are being widely studied and applied in chemical process engineering (Li et al. (2002), Li (2007), Li et al. (2008), Bartl (2006), Flemming et al. (2007), Henrion et al. (2001), Wendt et al. (2002)); optimal power flow (Shiina (1999)), finance and economics (Chanres and Cooper (1959)), water resource management (Prekopa (1995)), reliability based design optimization (Chan et al. (2006)), etc. Special cases (q2 = m) of (4) and (5) are either
(3)
where u, ξ, y are the vectors of the decision, uncertain input and uncertain output variables, respectively. The probability P r(·) is a measure of the confidence level of satisfying the inequality constraints; where α ∈ [0, 1] is a pre-given confidence (reliability) level. Stochastic optimization problems of this form are commonly known as chance (or probabilistic) constrained optimization problems. Assuming that f : Rn → R, g : Rn × Rm × Rl → Rq1 , and G : Rn × Rm × Rl → Rq2 , the chance constraint can have one of the following two forms: either P r{Gi (u, y, ξ) ≤ 0} ≥ αi , i = 1, . . . , q2 , αi ∈ [0, 1] (4) or P r{Gi (u, y, ξ) ≤ 0, i = 1, . . . , q2 } ≥ α, α ∈ [0, 1]. (5) The former are called separate or single chance constraints; while the latter are known as joint chance con⋆ We are very grateful to the DFG - Deutsche Forschungsgemeinschaft for supporting this research with Project No. LI806/8-2.
P r{yimin ≤ yi (u, ξ) ≤ yimax } ≥ αi , i = 1, . . . , m,
(6)
or P r{yimin ≤ yi (u, ξ) ≤ yimax , i = 1, . . . , m} ≥ α. (7) These constraints define the level of confidence that bounding constraints should be satisfied for a reliable as well as optimal operation of a system. The numerical computation of such chance constrained optimization problems is quite difficult and expensive. The usual approach is by relaxation of the problem into a deterministic non-linear optimization problem. 2. MONOTONICITY APPROACH TO THE SOLUTION OF CHANCE CONSTRAINED OPTIMIZATION PROBLEMS Of interest is the chance constrained nonlinear optimization problem (1) - (3) with either of the special chance constraints (6) or (7). It is assumed that the uncertain
variables ξ posses some joint probability density ρ(ξ). Such a density could be determined based on historical data of the application problem being considered. However, due to the (possible) nonlinear relationships between the variable y and ξ, it is usually hard to find out the probability density of the output variables y. Thus, the idea proposed by Wendt et al. (2002) is to transform the chance constraints in the space of y into chance constraints in the space of the variable ξ whose joint distribution is already known. The probability-measure preserving transformation can be guaranteed through a monotonic relation among the uncertain input variables ξ and output variables y. In case of separate chance constraints, it is assumed that for each of the restricted variables yi , there is a function F relating yi to one of the inputs ξj(i) such that ξj(i) 7→ fi (ξj(i) ) := F (ξ1 , . . . , ξj(i) , . . . , ξl , u) = yi
That is, there is a monotonic relation between yi and ξj(i) ; i.e., as ξj(i) increases, increases yi , too; ( it can also be the case that as ξj(i) increases, decreases yi ). There are several applied optimization models in which such monotonic relations readily available. (See the model in Section 4). In fact, the approach is to use the equality constraints g(u, y, ξ) = 0 in verifying such monotonicity relations between the variables yi and ξi ; furthermore, for a given u, to solve the systems of equations g(u, y max , ξ) = 0,
in order to determine ξ min (u) and ξ max (u), respectively; so that the (say, separate) chance constraints can be equivalently written as P r{yimin ≤ yi (u, ξ) ≤ yimax } min max = P r{ξj(i) (u) ≤ ξj(i) ≤ ξj(i) (u)} ≥ αi ,
for each i = 1, . . . , m. By doing so the equality constraints g(u, y, ξ) = 0 will be effectively eliminated from the optimization problem and the separate chance constrained problem will be transformed into min s.t. Z
f (u) +∞
...
Z
(8)
max ξj(i) (u)
...
min (u) ξj(i)
−∞
Z
+∞
ρ(ξ1 , . . . , ξl )dξl . . . dξ1 (9)
By defining hi (u) :=
+∞
... −∞
∇hi (uk ) =
Z
+∞
−∞ max k ∂ξj(i) (u )
...
Z
+∞
−∞
−
h
max k ρ(ξ1 , . . . , ξj(i) (u ), . . . ξl )×
min k ρ(ξ1 , . . . , ξj(i) (u ), . . . ξl )
min k ∂ξj(i) (u )
∂u
#
max ∂ξj(i) (uk ) ∂u
∂ξ min (uk )
and j(i) can be respecThe derivatives ∂u tively determined from the equation g(u, ξ, u) = 0 by the implicit function theorem. Consequently, the computation of h(uk ) and ∇hi (uk ) involves the evaluation of integrals. In fact, this approach requires efficient tools for numerical computation of multi-dimensional integrals. Properly speaking, the efficiency of the approach rests upon the precise computation of integrals. For this reason quadrature based methods, such as collocation on finite elements, are quite preferred instead of methods based on sampling. 3. DIRECTIONAL MONOTONICITY The above solution approach relies on a monotonic relationship between the input and output variables. Therefore, it is required to conduct a monotonicity analysis of a given system. However, such monotonic relations may not be available in all cases. Recently, monotonic properties of continuous dynamical systems are studied by Angeli and Sontag (2003). They give a general monotonicity property based on tangent cones; which they later on specialize by using partial derivatives. Here is given a similar, but weaker monotonicity property. Now, considering chance constraints (6) and for the sake of the practical problem to be considered in the next section, it is assumed that the constraint g(u, y, ξ) = 0 to have the form yj = gj (u, ξ1 , . . . , ξl ), j = 1, . . . , m; (i.e. m = q1 = q2 ).
−∞
≥ αi , i = 1, . . . , m, αi ∈ [0, 1].
Z
hi (u) ≥ 0, i = 1, . . . , m. This problem can be non-convex. Thus, an efficient numerical method is required for its solution. Currently, the method of choice for such nonlinear programming problem (NLP) is the interior-point method (see W¨achter and Biegler (2006))). However, to apply this NLP algorithm, at each step k, given uk , the constraint values hi (uk ) and gradients ∇hi (uk ) must be computed, where (see Xie et al. (2007), also Prekopa (1995))
dξl . . . dξj(i)−1 dξj(i)+1 . . . dξ1 .
ξj(i) = fi−1 (yi ).
and
f (u)
s.t.
∂u
and fi (·) is strictly monotonic; hence,
g(u, y min , ξ) = 0
min
Z
max ξj(i) (u)
...
min (u) ξj(i)
Z
+∞
ρ(ξ1 , . . . , ξl )dξl −∞
. . . dξ1 ) − αi , i = 1, . . . , m,
the problem (8)-(9) is written as
Definition 1. (Directional Monotonicity on a segment). Let g : Rn → R, y = g(ξ1 , ξ2 , . . . , ξn ) and a, b ∈ Rn , such that b > a . A function g is said to be monotonically increasing on the segment [a, b] in the direction of v = b−a iff the function ϕ(λ) = g(a + λv) is montonically increasing on [0, 1]. Theorem 2. Suppose g : Rn → R is a differentiable function and v ∈ Rn . Then g is monotonic increasing on the segment [a, b], a < b, if and only if ∇g(ξ)⊤ v ≥ 0, ∀ξ ∈ [a, b].
The proof follows trivially from the well known mean-value theorem (see for instance Spivak (1994)). The above directional monotony can be specialized to monotonicity in coordinate directions. In this case, let v be one of the standard unit vectors; i.e. v = ei (i ∈ {1, . . . , n}). Thus, ϕ(λ) = g(ξ + λei ) = F (ξ1 , . . . , ξi + λ, . . . , ξn ). As a result, g is monotonic on the segment [a, b] in the direction ei if and only if ϕ′ (λ) = λ∇g(ξ)⊤ ei ≥ 0, λ ∈ [0, 1]. Consequently, g is increasing on the segment [a, b] in the direction of the unit vector v = ei iff ∂g(ξ) ≥ 0, ∀ξ ∈ [a, b]. ∂ξi Accordingly, for y = g(ξ1 , ξ2 , . . . , ξn ), if it holds that ∂g(ξ) ∂ξi ≥ 0, then the variable y increases as the ith coordinate ξi increases. If y, ξ1 , . . . , ξn are random variables with some probability distributions, then it follows that P r {aj ≤ yj ≤ bj } = P r g −1 (ξ1 , . . . , ξi−1 , aj , ξi+1 , . . . , ξn ) ≤ ξi ≤ g −1 (ξ1 , . . . , ξi−1 , bj , ξi+1 , . . . , ξn )
The monotonicity properties being used in the next section can be interpreted as given above. Furthermore, with respect to the joint chance constraint (7), a monotonic relation should guarantee the transformation of a box constraint in the space of output variables y into a box constraint in the space of corresponding uncertain input variables ξ. In general, this is not a trivial issue. Considering the existing non-linear relation between y and ξ, work is still in progress to find out a viable monotonicity property.
Fig. 1. Buffer tank
a finite number of unit time steps t = 0, 1, . . . , N − 1, this problem is described as follows:
Objective f unction : φ(u0 , ..., uN −1 ) → min State equations : t = 1, . . . , N ; Vt = Vt−1 + qt−1 − ut−1 ; qt−1 Ct = Ct−1 + [COt−1 − Ct−1 ] ; Vt Initial conditions : V0 = v0 , C0 = c0 Control constraints : t = 1, . . . , N ; umin ≤ ut ≤ umax ;
4. AN OPTIMAL CONTROL PROBLEM WITH CHANCE CONSTRAINTS In a chemical process the feed to a distillation column can be an outflow of several upstream plants. For a smooth operation of the distillation column it is common to use a buffer tank to guarantee continuous supply to the column and avoid process disruption as shown in Fig. 1. Due to upstream processes there might be (uncertain) inflow disturbances to the buffer tank. The accumulation in the buffer tank should have certain minimum level (to guarantee continuity of the distillation process) and should not exceed the tank capacity. Moreover, the content of the buffer tank will be regulated by the outflow; i.e. the inflow to the distillation column. Let qt = q(t) represent the uncertain feed flow rate and COt = CO(t) inflow component concentration into the tank. The distillation column takes, at a time instant t, an amount ut = u(t) from the buffer tank for further processing. The volume Vt = V (t) and the concentration Ct = C(t) in the tank should satisfy bounding constraints up to a certain pre-given confidence levels α1 and α2 , respectively. The initial state of the tank, at t = 0, is given by V0 = V (0) and C0 = C(0). It is desired that the amount taken from the tank, from one time point to the next, must have a minimum fluctuation. Considering
Chance constraints : t = 1, . . . , N ; P {Vmin ≤ Vt ≤ Vmax } ≥ α1 ; P {Cmin ≤ Ct ≤ Cmax } ≥ α2 . Here, it is assumed that the uncertain inputs to have normal distributions qt ∼ N (µ(qt ), σt2 ) and COt ∼ N (µ(COt ), ̺2t ), where σ 2 and ̺2t are given variance of qt and COt , respectively, and µ(·) represents expected values. Linear problems of the above type have been considered, for instance, by Geibel and Wysotzki (2005), Henrion et al. (2001), Li (2007); Li et al. (2002); non-linear problems by Li et al. (2008) and Xie et al. (2007). In Li et al. (2002, 2008) a truncated normal distribution (see Horace (2005)) for the uncertain input variables has been used, while Henrion et al. (2001) studied the problem by using a scenario approach. In this paper, the uncertain variables are assumed to have standard normal distributions without truncation and integrals are evaluated by the direct numerical method. The approach used to solve the above problem is the moving horizon scheme. The control and prediction horizons are taken to be of equal length h. Thus, for each k = 0, . . . , N − 1, the following problem is to be solved:
min (∆u)⊤ ∆u s.t. i−1 X
56
(qk+j − uk+j ), i = 1, . . . , h
j=0
Ck+i
V − u k+j k+j Ck = V k+j+1 j=0 i−1 i−1 X Y Vk+s − uk+s qk+j COk+j , + V Vk+j+1 k+s+1 j=0 s=j+1
i−1 Y
54 C0(k) in g/l
Vk+i = Vk +
Different realizations of the concentration C0(k)
58
52 50 48 46 44
1
2
3
4
5
6
7
8
9 10 11 Time step k
12
13
14
15
16
17
18
19
i = 1, . . . , h;
umin ≤ uk+i ≤ umax , i = 1, . . . , h; P r{Vmin ≤ Vk+i ≤ Vmax } ≥ α1 , i = 1, . . . , h;
Fig. 2. Realization of the uncertain inflow concentration of COt
P r{Cmin ≤ Ck+i ≤ Cmax } ≥ α2 , i = 1, . . . , h;
=⇒ Vk+i ↑
and COk+i−1 ↑
=⇒ Ck+i ↑,
i = 1, . . . , h; meaning that Vt increases as qt−1 increases or Vt decreases if qt−1 decreases, etc. When transforming the chance constraint for Vk+i , the uncertain variables qk , . . . , qk+i−1 are taken into consideration; while w.r.t. Ck+i , the variables qk ,. . .,qk+i−1 , COk ,. . .,COk+i−1 are considered. Thus, in each horizon, it is assumed that each group of stochastic variables qk ,. . ., qk+i−1 and qk ,. . ., qk+i−1 , COk ,. . .,COk+i−1 have their respective multivariate normal distributions with corresponding covariance matrices. Accordingly, the length of the horizon determines the number and dimension of integrals to be computed. Furthermore, derivative values of constraints, required for the NLP algorithm, can be evaluated either iteratively using the method suggested by Prekopa for multivariate normal integrals (see Prekopa (1995) ) or by using quadrature methods (Bartl (2006)). In this paper a particular case of h = 3 for the purpose of numerical experimentation is considered. The following data are defined: • initial values V0 = 160, C0 = 50; • bounds Vmin = 130, Vmax = 170, Cmin = 49, Cmax = 51, umin = 0, umax = +∞; • reliability levels α1 = α2 = 0.8; • µ(qt ), µ(COt ), σt = V ar(qt ) = 0.7, ̺t = V ar(COt ) = 1, t = 1, . . . , N • N = 20. In each horizon the non-linear optimization problem is solved using IPOPT (W¨ achter and Biegler (2006)). Values of chance constraints and their gradients are computed using codes developed by Bartl (2006) based on orthogonal collocation. Integration limits are determined by solving the equality constraints using the Newton method. The following figures depict the results of numerical experiments for 21 realizations. Comparing Fig. 2 and Fig. 3, the realizations for COt and Ct justify the monotonic relations between them. It can
50.6 50.4 C(k) in g/l
qk+i−1 ↑
50.8
50.2 50 49.8 49.6 49.4 49.2 49
2
4
6
8
10 12 Time step k
14
16
18
20
Fig. 3. Resulting tank concentration for the 20 time intervals Different realizations of the feed q(k)
18 16 14 12 Feed q(k) in l
where ∆u = (uk+1 − uk , . . . , uk+h − uk+h−1 ), and define Qt−1 s=t (·) = 1. It is easy to observe that the following monotonic relationships hold true:
Resulting tank concentration C(k) for realizations of q(k) and C0(k) and calculated u(k)
51
10 8 6 4 2 0
1
2
3
4
5
6
7
8
9 10 11 Time step k
12
13
14
15
16
17
18
19
Fig. 4. Realization of the uncertain feed flow qt also be seen from Fig. 4 and Fig. 5 that the volume Vt of the buffer tank relates with the inflow qt monotonically. It is not surprising to observe from Fig. 5 and Fig. 6 that the control ut imitating the behavior of the tank volume Vt : when there is little in the tank, little will be allowed to flow out; and when there is more, then more will be allowed to flow out. The realized values of tank volume Vt (Fig. 5) and concentration Ct (Fig. 3) satisfy, in almost all the 21 realizations, their bounding constraints up to the required reliability levels α1 = 0.8 and α2 = 0.8. Good choice for α1 and α2 are values greater or equal to 0.5 (at least 50% reliability). Surely a process is expected to work with a higher relia-
REFERENCES
Resulting tank volume V(k) for realized q(k) and calculated u(k)
180
170
V(k) in l
160
150
140
130
120
2
4
6
8
10 12 Time step k
14
16
18
20
Fig. 5. Realized tank volume Resulting controls u(k) for the realizations of q(k) and C0(k)
18 16 14
u(k) in l
12 10 8 6 4 2
1
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16 17 18 19 20 Time step k
Fig. 6. Control realization for the 20 time intervals bility level; i.e. with values for α1 and α2 near 1. However, requiring higher reliability levels may cause infeasibility of the constraints. The analysis of the largest values of α1 and α2 that still guarantee the feasibility of the problem has not been done yet. This will be a future undertaking. A similar investigation for linear systems can be found in Li et al. (2002). The numerical experiments in this work were carried out using - Matlabr R2008a for Linux 64-bit, - IPOpt Version 3.4.2 with the linear solver ma27, - Processor: Athlon64 X2 3800+, 2GB RAM. Integrals are calculated with a maximum absolute error of 0.002 and the NLP algorithm is allowed to work with a termination tolerance of 0.0001. Individual realizations required an average total computing time of approximately 100 minutes. Among these, most of the time spent went for the calculation of integrals. CONCLUSIONS: Nonlinear dynamic optimization under chance constraints is considered. The solution of the problem depends on a monotonicity relation between uncertain input and output variables. For this purpose a monotonicity analysis is introduced to transform chance constraints from the space of uncertain output variables into the space of uncertain input variables. The moving horizon scheme is used for the optimal control of a tank process under the feed flow and concentration uncertainties.
Angeli, D. and Sontag, D. (2003). Montone control systems. IEEE Trans. on Auto. Ctrl., 48, 1684–1698. Bartl, M. (2006). Station¨are Prozessoptimierung unter Einbezug von unsicheren Randbedingungen am Beispiel einer Destillationsanlage. Diplomarbeit, Fakult¨ at f¨ ur Informatik und Automatisierung, Technische Unviersit¨ at Ilmenau. Chan, K.Y., Skerlos, S.J., and Papalambros, P.Y. (2006). Optimal design with non-normally distributed parameters, conditional probability, and joint constrained reliabilities. Proceedings of the 2006 ASME International Design Engineering Technical Conferences and Computers and Information In Engineering Conference, Philadelphia, PA, DETC2006/DAC-99417, 14 Pages. Chanres, A. and Cooper, W. (1959). Chance-constrained programming. Management Science, 6, 73–79. Flemming, T., Bartl, M., and Li, P. (2007). Set-point optimization for closed-loop control systems under uncertainty. Industrial & Engineering Chemistry Research, 46, 4930–4942. Geibel, P. and Wysotzki, F. (2005). Risk-sensetive reinforcement learning apllied to control under constraints. J. Artificial Inteligence Res., 24, 81–108. Henrion, R., Li, P., M¨oller, A., Seinbach, M.C., Wendt, M., and Wonzy, G. (2001). Stochastic optimization for operating chemical process under uncertainity. In M. Gr¨othschel, S. Krumke, and J. Rambau (eds.), Online Optimization of Large Scale Systems, 455 – 476. Springer Verlag. Horace, W.C. (2005). Some results on the multivariate truncated normal distribution. Journal of Multivariate Analysis, 94, 209–221. Li, P. (2007). Prozessoptimierung unter Unsicherheiten. Oldenbourg Verlag. Li, P., Arellano-Garcia, H., and Wozny, H. (2008). Chance constrained programming approach to process optimization under uncertainity. Computers and Chemical Engineering, 32, 24–45. Li, P., Wendt, M., Arellano-Garcia, H., and Wozny, G. (2002). Optimal operation of distillation process under uncertain inflows accumulated in a feed tank. AIChe Journal, 48, 1198–1211. Prekopa, A. (1995). Stochastic programming. Kluwer Academic Publishers. Shiina, T. (1999). Numerical solution technique for joint chance-constrained programming problem- application to electric power capacity expansion. J. of OR Socity of Japan, 42, 128 – 140. Spivak, M. (1994). Calculus, 3rd ed.,. Publish or Perish Inc. W¨achter, A. and Biegler, L.T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Prog., 106, 25 – 57. Wendt, M., Li, P., and Wozny, G. (2002). Nonlinear chance-constrained process optimization under uncertainity. Ind. Eng. Chem. Res., 41, 3621 – 3629. Xie, L., Li, P., and Wozny, G. (2007). Chance constrained nonlinear model predictive control. In R. Findeisen, F. Allg¨ower, and L.T. Biegler (eds.), Assessment and Future Directions of nonlinear model predictive control, 295 – 304. Springer Verlag, London.