Journal of Process Control 12 (2002) 373–389 www.elsevier.com/locate/jprocont
Real-time optimization under parametric uncertainty: a probability constrained approach Yale Zhang, Dayadeep Monder, J. Fraser Forbes* Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB, Canada, T6G 2G6 Received 6 September 2000; received in revised form 4 June 2001; accepted 24 August 2001
Abstract Uncertainty is an inherent characteristic in most industrial processes, and a variety of approaches including sensitivity analysis, robust optimization and stochastic programming have been proposed to deal with such uncertainty. Uncertainty in a steady state nonlinear real-time optimization (RTO) system and particularly making robust decisions under uncertainty in real-time has received little attention. This paper discusses various sources of uncertainty within such closed loop RTO systems and a method, based on stochastic programming, that explicitly incorporates uncertainty into the RTO problem is presented. The proposed method is limited to situations where uncertain parameters enter the constraints nonlinearly and uncertain economics enter the objective function linearly. Our approach is shown to significantly improve the probability of a feasible solution in comparison to more conventional RTO techniques. A gasoline blending example is used to demonstrate the proposed robust RTO approach. # 2002 Elsevier Science Ltd. All rights reserved. Keywords: Real-time optimization; Parametric uncertainty; Stochastic programming; Joint probability constraints
1. Introduction Real-time optimization (RTO), or on-line operations optimization, has a wide appeal in the process industries because of its promise for improving process profitability. Unfortunately, many of these systems have poor long-term service factors due, at least partially, to inadequate optimization robustness [1]. RTO is a model-based process control approach that uses current process information (i.e. process model and economic data) to predict the optimal operating policy for a process unit during the next RTO interval. Although most RTO systems attempt to improve model accuracy through model updating, there is a large number of uncertainty sources in the RTO problem and these uncertainties can seriously degrade the performance of any RTO system. A key issue in RTO is ensuring that any predictions will result in feasible operation given the recognized uncertainty in the RTO problem. In this work the term robustness of a model-based RTO system speaks to the * Corresponding author. Tel.: +1-780-492-0873; fax: +1-780-4922881. E-mail address:
[email protected] (J. Fraser Forbes).
ability of the system to predict an ‘‘optimal’’ operating policy that will be feasible at some level of probability. Then, a robust RTO system will usually provide a more conservative solution than a conventional RTO system to ensure that process operating constraints are unlikely to be violated. Such a robust RTO system ‘‘trades-off’’ some potential profitability for an assurance of the feasibility of the predictions. Uncertainty is inherent in most process control and optimization problems. It is well recognized that uncertainty in the problem data (e.g. noisy, incomplete, unmeasurable or erroneous data) is an essential consideration in both process design/synthesis and process operation optimization. Such uncertainty was classified by Kraslawski [2] into two types: (1) ambiguity, if its value cannot be definitely established; and (2) imprecision, if its value is not sufficiently determined with respect to a given scale. In that work, the author discussed the effect of various uncertainties in process modeling, control and synthesis. Uncertainty in process design/planning was also investigated by Ierapetritou et al. [3]. A number of approaches have been proposed to deal with the uncertainties in engineering optimization applications, including: (1) solution sensitivity/stability
0959-1524/02/$ - see front matter # 2002 Elsevier Science Ltd. All rights reserved. PII: S0959-1524(01)00047-6
374
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
Nomenclature b D d e f g H h i,j K k J m N n R r u x z l F
=constant vector =deterministic set =constant vector =error between setpoint and process measurement =equality constraints (process models) =inequality constraints =parameter matrix in linearized models =parameter vector in linearized models =index =feasible domain =number of constraints =objective function =number of dependent variables =Normal distribution =number of decision variables =correlation coefficient matrix =controller setpoint =dependent variables, or controller output =decision variables =process measurements =constant =adjustable model parameters (model bias) =expected value =standard deviation =random variable =uncertain economic information in objective function =parametric uncertainty in constraints =probability level =constant =constant =probability distribution function
analysis; (2) robust optimization, and (3) stochastic optimization. Sensitivity analysis was originally employed to measure the impact of changes in the problem data (i.e. parametric uncertainties) on the optimal solution. In order to achieve solutions that are less sensitive to variations in the problem data, Becker et al. [4] developed a robust optimization approach by augmenting the objective function of the optimization problem with a penalty term intended to minimize parametric sensitivity. Robust optimization [5,6] is used to optimize the expected value of a chosen performance index for a given level of risk and includes formulations based on expectation, weighted mean-variance, worst-case and so forth. Alternatively, stochastic optimization [7,8] attempts to directly solve a problem given uncertainty in the data. It handles these uncertainties using recourse functions, chance constraints and so forth, to transform the stochastic optimization problem to an equivalent
R
=real numbers
Operators E[.] =expectation operator Var[.] =variance P[.] =probability r =gradient Subscripts =original RTO system or current operating 0 point in =inner out =outer Superscripts =iteration index * =optimum point ˆ =expected value ˆ =predicted or estimated
(i)
Gasoline blending operation study A,B =constant matrix be,bu =constant vector cf =cost of feedstock cp =cost of product e =vector containing ones (i.e. e=[1,...,1]T) qf =feedstock quality s =product specifications x =feedstock flow rate for regular grade gasoline blending y =feedstock flow rate for premium grade gasoline blending t =time C =constant matrix " =noise vector
deterministic optimization problem, which can be solved using a number of readily available optimization techniques. Since stochastic optimization deals directly with uncertainty in an explicit, ‘‘natural’’ manner, it is a good candidate technology for use in RTO. Although there has been considerable effort in the areas of off-line optimization under uncertainty, making robust decisions under uncertainty in real-time has received little attention. Consider an RTO system that includes components for detecting steady-state operation, gathering and validating plant data, updating the process model, estimating optimum plant operation, analyzing optimization results, and transferring optimal operating policies to advanced controllers, there are often many sources of uncertainty around this closed RTO loop, which influence the model-based RTO predictions. This paper begins with discussing four sources of uncertainties (i.e. process uncertainty, market uncertainty,
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
measurement uncertainty and model uncertainty) in the RTO system and their influences on RTO predictions. This is followed by a brief overview of the existing approaches to deal with these uncertainties. Then, a robust RTO formulation is presented that is able to account for various uncertainties in the closed RTO loop, to make robust decisions under such uncertainty in real-time, as well as to increase the operating profit of the plant. The proposed approach is limited to problems where parameters enter the constraints nonlinearly and the economic uncertainty enters the objective function linearly. This formulation addresses a substantial fraction of the RTO problems encountered in industry. Although the main focus of this paper is the development of a practical formulation for the Robust RTO problem and a solution method for this formulation, other contributions include a clear discussion of how the process operating constraints are naturally cast as probabilistic constraints and how these constraints can either be treated independently or jointly within the optimization problem. Further, the benefits and effects of treating the probabilistic constraints as independent or correlated are compared. The proposed approach is illustrated by a gasoline blending example that is taken from Forbes and Marlin [10] and modified slightly for the purposes of this paper.
2. RTO under uncertainty A typical model-based RTO loop is shown in Fig. 1. First, the plant data (z) are gathered, validated to avoid gross errors in the process measurements, and may be reconciled using material and energy balances to ensure
375
the data set used for model updating is self-consistent. These validated measurements (z0 ) are used to estimate some set of adjustable model parameters () to ensure the model represents the plant, as accurately as possible, at the current operating point. Then, the optimum controller setpoints (xm ) are calculated using the updated model, and are transferred to the advanced process controllers after they are checked by the command conditioning subsystem. It is worth noting that the RTO system in this work is employed to determine optimal steady state operating points in real-time, using rigorous, nonlinear process models and available economic information so as to maximize process profitability. This paper does not consider the on-line dynamic optimization problem contained in nonlinear model predictive controllers. Further, in this paper it is assumed that the underlying control system will implement the calculated setpoints, while ensuring a safe and feasible operation. An RTO problem is subject to four types of uncertainty [11]: 1. Market uncertainty, which refers to imprecise knowledge of the process economics, raw material availabilities, product demands, etc., 2. Process uncertainty, which refers to imprecise knowledge of the operation due to various process disturbances including unmeasured variations in ambient conditions, feedstock qualities, etc., 3. Measurement uncertainty, which refers to imprecise knowledge of measured process variables due to sensor errors consisting of random measurement noise, biases, etc., 4. Model uncertainty, which refers to imprecise model predictions due to plant/model structural and parametric mismatch.
Fig. 1. RTO loop with uncertainties.
376
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
The variations that are the underlying cause of both process and market uncertainties may be either stationary or nonstationary and are often considered to be stochastic in nature. These two uncertainty sources affect both the objective function and constraints in the process optimization problem. Therefore, these two uncertainty sources can result in an RTO solution that is infeasible with respect to process operation, if they are not dealt with appropriately. Measurement uncertainty is the result of gross errors and random measurement noise. Gross errors includes process leaks, instrumentation biases, and so forth, which may happen unpredictably. Random measurement noise is often considered to be normally distributed with zero mean and a given covariance structure, and to have no correlation in time. In general, the latter occurs at high frequency (compared to the process changes to which the RTO system is expected to respond) during the normal operation of the plant, and may not influence the optimal operations; however, such high frequency variation in the process measurements may result in unnecessary, frequent changes in RTO predictions due to transmission of the random variation around the closed RTO loop (see Fig. 1). The Command Conditioning subsystem should eliminate transmission of such high frequency noise and the resulting variation in the RTO results, thereby increasing the plant operating profit. Plant/model mismatch is a crucial issue to any successful model-based RTO system. Structural mismatch arises from approximation of the underlying process mechanisms; whereas, parametric mismatch may stem from some process changes such as exchanger fouling, catalysts deactivation, etc. Then, model uncertainty may lead to RTO solutions that are at a significant distance from the true plant optimum, or even an infeasible solution with respect to the actual process constraints. In most RTO systems, model uncertainty is addressed through some form of on-line model updating. The three methods that are used to deal with these uncertainties in RTO to date include: (1) the back-off approach; (2) results analysis; and (3) stochastic optimization. In the back-off approach, developed by Loeblein et al. [12], a back-off vector is introduced into the optimization constraints to help ensure that the RTO system produces feasible operating points. Both model and measurement uncertainty in the closed RTO loop are taken into account to determine the value of the back-off term, which is calculated off-line and is not part of the RTO calculations. The main disadvantage of this approach is the inaccuracy caused by infrequent calculation of the back-off term, which is expected to change with operating level and active constraint sets. Results analysis in RTO systems, developed by Miletic and Marlin [13] and later extended by Zhang et al.
[11], uses multivariable statistical hypothesis testing to decide whether RTO predictions represent statistically significant changes with respect to the current plant operation, or are the result of the transmission of random measurement noise around the RTO loop. The results analysis methods have been shown to be effective in reducing the influence of measurement noise on RTO predictions. Unfortunately, results analysis methods do not explicitly address the effect of measurement uncertainty on the RTO system, as these approaches are postoptimization calculations. Of the available approaches, stochastic optimization is the only strategy that explicitly incorporates uncertainty into the solution of the operations optimization problem. It does so by using recourse functions, chance constraints and so forth, to transform the stochastic program to an equivalent deterministic optimization problem. Monder and Forbes [14] used chance constrained programming in the development of a robust time horizon based RTO system for gasoline blending process to increase the operating profit. However, the chance constraints in [14] were treated separately (i.e. the covariance of uncertainties or correlation among the process constraints was ignored). Schwarm and Nickolaou [15] incorporated chance constraints into a linear model predictive controller (MPC) to deal with some forms of parametric uncertainty. The resulting controller improved the robustness of MPC with respect to output constraint satisfaction. Later Li et al. [16] extended this work by considering a robust MPC problem under joint probability constraints. In both of their approaches, the economic performance and market uncertainty were not addressed. This issue was discussed by Kassmann et al. [17] in their recent work. A two-tier MPC algorithm, a steady-state target calculation followed by a dynamic optimization, was developed. In their method, the target calculation was formulated as a linear stochastic programming problem to account for parametric uncertainty in the steady-state gain matrix of the nominal operation point; however, such steady-state gain models may not represent the process correctly when the nominal operation point is significantly changed. Although the MPC controller is still stabilizing for large changes in the operation, the results of the economic optimization may no longer be either optimal or feasible. Thus, as a number of researchers have indicated, stochastic optimization seems to deal with the different uncertainties within the closed RTO loop in the most comprehensive way. It provides approaches for dealing with parametric uncertainties in both objective function and constraints, and the uncertain constraints can be handled independently to simplify solution procedure, or handled jointly to address any correlation among the constraints. In the next section, a robust RTO formulation is proposed and solved using a stochastic optimization approach.
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
3. Robust RTO Most RTO problems can be cast as: min
ðx;uÞ2D0
J0 ðx; u; Þ
s:t: : f0 ðx; u; Þ ¼ 0 g0 ðx; u; Þ 5 0
ð1Þ
where: x Rn is a vector of decision variables; u Rm is a vector of dependent variables; H is a vector representing the uncertain economic information (market uncertainty) in the objective function; U is a vector representing the parametric uncertainty (process uncertainty, parametric model uncertainty, measurement uncertainty, etc.) in the process constraints; J0 is a nonlinear objective function and is often the operating profit in RTO; f0 and g0 are sets of equality constraints (i.e. the process model) and operating constraints that contain parametric uncertainties, respectively, and are nonlinear in x,u and ; and D0 is a deterministic set usually defined by a finite number of inequalities involving functions of x and u. In RTO applications where the dependent variables can be eliminated using the process models, an equivalent RTO problem in the reduced space of the decision variables can be represented as: min x2D
Jðx; ; Þ
s:t: : gðx; Þ 5 0
ð2Þ
Note that in problem (2) the subscript ‘‘0’’ has been eliminated to indicate the reduced space representation of problem (1). The developments in this work are limited to RTO formulations that can be represented as in problem (2).
objective function J with respect to the various uncertainties (i.e. and ); and the second term Var½Jðx; ; Þ
represents the ‘‘spread’’ of the distribution of objective function values about the expected value. The solution of problem (3) results in a trade-off between the expected value and the variance of the objective function (i.e. an operating point that is the best compromise between profitability and expected variation of that profit). The choice of the constant is arbitrary and reflects the decision maker’s opinion as to the relative importance of the expected value and the variance of the optimum objective function value. The variation in RTO parameters is often assumed to be reasonably approximated by a normal distribution and as a result, can be completely characterized by the first two moments of the distribution. Further, RTO systems are closed-loop, steady-state controllers that attempt to locate the optimum operation though a series of small moves (i.e. usually move limits are placed on the RTO outputs). Then, the variance term in problem (3) may be considered approximately constant from one RTO interval to the next and the problem reduces to: min E½Jðx; ; Þ
x2D
s:t : gðx; Þ 5 0
ð4Þ
The solution of problem (4) results in a set of decision variable values which produce the minimum objective function value, on ‘‘average’’. Generally, economics enter into the objective function linearly (i.e. @J @ is independent of ); however, model parameters do not. Since the value of the uncertain model parameters can usually be assumed to be independent of the uncertainties in economics and at any RTO interval, is assumed ^ proto be constant at the current estimated value vided by on-line model updating, problem (4) becomes: ^ ; Þ
min E½Jðx; x2D
s:t: :
3.1. Objective function uncertainty
gðx; Þ 5 0 There are several approaches for dealing with uncertainty in the objective of an optimization problem. One possible formulation [5,6] can be taken as a linear combination of the expectation and the standard deviation of the objective function Jðx; ; Þ: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi min E½Jðx; ; Þ þ Var½Jðx; ; Þ
s:t: : gðx; Þ 5 0
377
ð3Þ
where 04a41 is a user specified constant. The first term E½Jðx; ; Þ represents the expected value of the
ð5Þ
Futher, since the objective function of RTO problem is assumed to be linear in the economic parameters (), then problem (5) reduces to: ^ ; Þ min Jðx; x2D
s:t: : gðx; Þ 5 0
ð6Þ
where: is the expected value of the economic parameters. Thus, in most RTO applications it may be sufficient to use the expected values of the problem
378
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
parameters to deal with the effect of uncertainty in the objective function. For the purposes of this work, the objective function will be evaluated in terms of the expected values of the problem parameters.
Pðgj ðx; Þ 5 0Þ 5 j ðj ¼ 1; ; kÞ
has received much more attention in the literature. Thus, the robust RTO problem can be written as either an individual probability constrained (IPC) problem:
3.2. Uncertain constraints min There are two approaches to formulating an optimization problem when there is uncertainty in the constraints. These are expectation-based constraints and probability constraints [8]. Expectation constraints can take either the simple form: E½gðx; Þ 5 0
ð7Þ
x2D
^ ; Þ Jðx;
s:t: : Pðgj ðx; Þ 5 0Þ 5 j ðj ¼ 1; ; kÞ
min
^ ; Þ Jðx;
s:t :
or the conditional form:
Pðgðx; Þ 5 0Þ 5 0
ð12Þ
ð8Þ
where d is a limitation for the average measure of constraint violation, if the uncertain constraints are violated (i.e. gðx; Þ < 0Þ. In the first case, the constraints are expected to be satisfied ‘‘on average’’ but the amount of constraint violation is not restricted, and the constraints can be expected to be violated no more than approximately 50% of the time. In order to limit the amount of constraint violation, Inequality (8) requires that whenever the constraints are violated, the ‘‘average’’ violation must be within a specified bound d. In bounding the ‘‘average’’ constraint violation, inequality (8) ignores the frequency of constraint violation. Then, in order to control the magnitude of violation and gain some control of the frequency of violation, both constraint types would be required. In the many industrial applications, a more ‘‘natural’’ manner framework for discussing uncertainty in the constraints may be to pose the problem in terms of the probability of constraint violation, which is the approach taken in this work. For a given probability distribution, expression of the constraints in this manner would provide some control over both the frequency and magnitude of the variation. Then, within a probabilistic framework, a well-designed RTO system would ensure that the constraints are respected to some acceptable level of probability. This paradigm would yield a probability (chance) constrained [8,18] RTO problem. A probability constraint can be stated as: Pðgðx; Þ 5 0Þ 5 0
ð11Þ
or a joint probability constrained (JPC) problem: x2D
E½ gðx; Þjgðx; Þ < 0 4 d
ð10Þ
ð9Þ
Inequality (9) states that all of the process constraints have to be satisfied at or beyond a specific level of probability (0), simultaneously. Such joint probability constrained problems have proven very difficult to solve and as a result, the individually constrained problem:
Although both problems (11) and (12) attempt to achieve an optimal solution, which ensures that the inequality constraints are satisfied at some level of probability, there is a distinct difference between them. In the JPC problem, satisfaction of the whole set of uncertain constraints is represented by one probability constraint with the risk level of 0, and it is guaranteed that the solution is feasible with respect to all of the constraints and the uncertainties therein at not less than the 0 probability level. Whereas, each individual constraint in the IPC problem is satisfied at a certain probability level (i); however, the probability level at which all of the constraints are satisfied is less than any one of individual probability levels. For example, consider a certain solution of an IPC problem with 10 probability constraints. Each probability constraint is assumed to be satisfied at or beyond the probability level of 95% (i.e. j ¼ 95%; j ¼ 1; ; 10). Then, the worst case (i.e. lower bound) probability level at which this solution satisfies all of the constraints is equal to: 10
0 ¼ j ¼ ð95%Þ10 60%. A comparison of the IPC j¼1
and JPC via simulation studies will be given in Sections 4 and 5. A crucial issue that must be considered when employing probabilistic constraints is ensuring that the constraint satisfaction probability (i) is not set so high that no feasible solution exists for the problem. (Consider that as i!100% for a specific probabilistic constraint, the required ‘‘safety margin’’ to ensure that gi ðx; Þ 5 0 is respected goes to 1). The rate at which the size of the ‘‘safety margin’’ grows with depends on the variance term in the uncertainty description for the parameters () appearing in the constraint. Thus, the constraint satisfaction probabilities cannot be arbitrarily specified and care must to taken when setting these probability levels. Finally, should the situation be encountered where there is no feasible solution at a
379
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
specified constraint satisfaction level, then either the constraint satisfaction probability must be relaxed (i.e. reduce i) or effort must be expended to reduce the amount of uncertainty in the appropriate model parameters.
observation ðx0 ; 0 Þ, the probability constraints can be rewritten as: PðhTj x 5 j Þ 5 j
ð13Þ
where: 4. Solution method As previously developed, a robust RTO problem can be cast as either IPC or JPC problem when the probability constraints are employed. It is worth noting that the inequality constraints represented by gðx; Þ 5 0, which are contained within the probabilistic constraints in problems (11) and (12), may be either linear or nonlinear. Currently, there is very little stochastic optimization literature for situations where the uncertain parameters enter nonlinearly into the problem statement. The method proposed in this work linearizes the constraint set and uses ‘‘bias’’ updating to compensate for plant/model mismatch. It has been noted that although the bias update technique is widely used in process industries for the purposes of process control and optimization, this form of model updating cannot guarantee convergence to the true plant optimum [9]. It is our contention that in an uncertain RTO problem, it is more important to find a solution that is expected to satisfy all of the operating constraints, given the specified uncertainties, than it is to calculate optimal operations based on the nominal values of the problem parameters. In other words, if the uncertainties are significant in either the objective function or constraints of an optimization problem, then the feasibility of optimal solution should be the primary concern in solving such a problem. Thus, this work is applicable to process models that pass the model accuracy test given by Forbes and Marlin [9], which ensures that model updating technique can provide the necessary values of ‘‘bias’’ terms from available process measurements such that the true plant optimum can be predicted, and to RTO problems where a high probability of solution feasibility is of primary importance. Finally, the linearization approach that is proposed here may result in poor or infeasible decisions, when substantial nonlinearity is present in the plant. To ensure that a solution is feasible, a posterior simulation using the nonlinear model through sampling the uncertain parameters may be used to check the range of the discrepancy between the nonlinear and linearized model. Alternatively, infeasibility could be caught in the command conditioning (or results analysis) component of the RTO system or by the underlying control system. 4.1. IPC algorithm By employing a first-order Taylor series approximation for gj ðx; Þ in the IPC problem (11) at the current
hj ¼ ðrx gj Þx0 ; 0 j ¼ r T gj
x0 ; 0
ð14Þ þ j þ bj
ð15Þ
and j is a bias term that is used to compensate for the plant/model mismatch including the effects of model linearization. It can be considered as a random variable, since it is updated using uncertain process measurements in every RTO interval. Further, bj is a constant given by: T gj 0 gj ðx0 ; 0 Þ ð16Þ bj ¼ ðrxT gj Þx0 ;0 x0 þ r x0 0
When both and follow the multivariate normal distribution, then: j Nð ; 2j Þ
ð17Þ
where the mean value ( j ) and variance ( 2j ) can be calculated from the statistical properties of , and the constant bj. Assuming that the probability distribution function of j is given by Fj , then as in [7]: PðhTj x 5 j Þ ¼ Fj ðhTj xÞ 5 j
ð18Þ
that is: hTj x 5 F 1 j ðj Þ
ð19Þ
Hence, the IPC problem (11) is equivalent to: min x2D
^ ; Þ Jðx;
s:t: : hTj x 5 F 1 j ðj Þ ðj ¼ 1; ; kÞ
ð20Þ
and can be solved using a wide range of available optimization algorithms. 4.2. JPC algorithm Prekopa and Szantai [19,20] presented an algorithm based on the supporting hyperplane method [21]. This is the approach adopted for solving the JPC problem in this paper. Again, as in the IPC algorithm, the nonlinear constraints are linearized and bias updating is employed, and the resulting JPC problem is:
380
min x2D
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
^ ; Þ Jðx;
^ ; Þ min Jðx; x2D
s:t: :
s:t: : PðHx 5 Þ 5 0
x 2 Kðiþ1Þ
ð21Þ
ðiþ1Þ is an If PðHxðiþ1Þ out 5 Þ 5 0 , then xout optimal solution to problem (21); otherwise, go to step (a), using (i+1) instead of (i).
where H ¼ ½h1 h2 hk T and ¼ ½ 1 2 k T . Then, the JPC algorithm consists of two phases: Phase 1. Find two initial solutions. One is feasible, that is, xin 2 Kð1Þ ¼ xjP½Hx 5 5 0 ; and the ð1Þ other is an infeasible solution, xð1Þ out 62 K . To obtain xin, the following linear programming problem is solved [20]: min x2D
k X
hi1 x1 þ hi2 x2 þ þ hin xn i = i
ð26Þ
The algorithm can be speeded up if a fixed xin in Phase 2 is replaced by an adaptive term. Pre´kopa [8] showed that the choice: xðiþ1Þ ¼ xðiÞ in in þ
i¼1
1 ðiÞ xl xðiÞ in iþ1
ð27Þ
s:t: : Hx 5 i þ i ð22Þ where the recommended value for is 3, as given in the literature. The strategy of adjusting the value of is also given in [20]: If the optimal solution to problem (22) is out of the feasible domain K(1), then the value of should be enlarged; however, a relatively large value may make the problem (22) too restricted to solve. An effective way to obtain xð1Þ out is to solve the corresponding IPC problem (20), as a feasible solution to an IPC problem often does not satisfy the joint probability constraint. Phase 2. There are two steps in phase 2 for iteration i: Step (a) Let l(i) be the largest l(04l41), for which x ð23Þ P H xin þ l xðiÞ 5 5 0 in out
are efficient in many cases. Note that in the supporting hyperplane method it is necessary to calculate the probability values and gradient vector of the probability as a function of the variables x at the actual boundary point xðiÞ l . For a multivariate normal distribution, given the mean values and covariance matrix of the random variables, a Monte Carlo approach is often used for probability calculations and is adopted for the purposes of this work. Care must be exercised in using the Monte Carlo approach to ensure that a sufficient number of simulations are performed to ensure an accurate approximation of the underlying distribution is obtained for optimization purposes. 4.3. Illustrative example A simple analytical example taken from [8] is selected in this paper to illustrate the IPC algorithm, JPC algorithm and their differences. The stochastic optimization problem is as follows: min x1 þ 2x2
is satisfied. Denote ðiÞ xðiÞ ¼ x þ l x x ð24Þ in in out l ðiÞ If PðHxðiÞ l 5 Þ P Hxout 5 < ", where " is a user specified tolerance, then xðiÞ l is an approximate solution to problem (21). Otherwise, define Kðiþ1Þ ¼ n o Þ50 xjx 2 KðiÞ and ðrx PðHx 5 ÞÞxðiÞ ðx xðiÞ l l
ð25Þ
x2D
s:t: : x1 þ x2 5 1 þ 2 x2 5 2
½1
½2
ð28Þ
where D ¼ fxj0 4 x1 4 2; 0 4 x2 4 5g; 1 and 2 follow a joint normal distribution with expectation m=[1,2]T, standard deviation s=[0.1, 0.2]T, and correlation coefficient matrix R¼
1 0:8
0:8 1
ð29Þ
then go to Step (b) Step (b) Let xðiþ1Þ be an optimal solution to the folout lowing problem:
In the simulations, the acceptable probability level for the satisfation of the uncertain constraints in either IPC or JPC problems is taken as 90%:
381
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
x1 þ x2 5 1 þ 2 0 ¼ P 5 0:90 x2 5 2
ð30Þ
and 1 ¼ Pðx1 þ x2 5 1 þ 2 Þ 5 0:90
ð31Þ
2 ¼ Pðx2 5 2 Þ 5 0:90
ð32Þ
The optimal solutions are given in Table 1. Table 1 shows that the deterministic solution is not robust to the uncertainties in 1 and 2 . The probability of satisfying both constraints 1 and 2 is 42.23%, and as expected, the probabilities of satisfying each individual constraint are approximately 50% each. Clearly, the IPC approach improves solution robustness considerably and as specified, provides a probability of satisfying each of the individual constraints at a level of approximately 90%; however, the probability of satisfying both constraints 1 and 2 jointly is only approximately 87%. This is due to the correlation between the constraints and the inability of the IPC algorithm to handle this correlation. In this example, the joint constraint satisfaction probability delivered by the IPC approach is only marginally less than the specified level of 90%. This is due to the comparatively small size of the problem. In general, the joint constraint satisfaction probability level obtained using the IPC algorithm will depend on two main factors: the number of probability constraints and the amount of correlation among them. The JPC algorithm delivered the specified 90% constraint satisfaction probability at the expense of inflating the objective function to 5.4804, which is only a little larger than the value of 5.4282 delivered by the IPC algorithm and the value of 5 obtained by treating the problem as deterministic at the expected values of the model parameters. As previously mentioned, the number of Monte Carlo simulations used in the solution of the JPC problem plays a significant role in the accuracy of the optimization results. To illustrate this, the JPC form of problem (28) was solved using three different numbers of Monte Carlo simulations: 500, 5000 and 50000. Then, the JPC Table 1 Optimal solutions for illustrative example
Optimization variable x1 Optimization variable x2 1 2 0 Objective function value
Deterministic solutiona
IPC solution
JPC solution
1 2 49.79% 49.74% 42.23% 5
0.9156 2.2563 90.12% 90.04% 86.88% 5.4282
0.9213 2.2795 93.29% 91.90% 90.01% 5.4804
form of problem (28) was solved 50 times using different sample realizations for each number of Monte Carlo simulations. These simulation results are shown in Fig. 2. As expected, the variance of the JPC solution decreases as the number of Monte Carlo simulations used in the optimization increases. Fig. 2 also shows the differences between the JPC solutions and the IPC solution.
5. Gasoline blending case study The example considered here is a gasoline blending process that is taken from [9] and modified slightly for the purpose of this paper. The process, shown in Fig. 3, involves simultaneously blending two grades (i.e. regular and premium) of automotive gasoline using five feedstocks: LSR naphtha, reformate, n-butane, catalytic gas, and alkylate. The optimization objective is to determine the feedstock flow rates at each RTO interval that maximize the profit, while satisfying the product specifications for each gasoline grade. These product specifications include: (1) minimum limits on research octane number (RON) and motor octane number (MON); and (2) a maximum limit on Reid vapor pressure (RVP). Other operating constraints include bounds on product demands and feedstock availabilities. The RTO problem can be written as: T e x max cTp T cTf ðx þ yÞ x;y e y ð33Þ ðaÞ s:t: : fðx; y; q f Þ 5 s bl 4 Ax þ By 4 bu
ðbÞ
where: x; y 2 R5 are decision variables representing the feedback flow rates for regular and premium grades gasoline blending, respectively; cp=[33, 37], cf=[34, 26, 10.3, 31.3, 37] are economic values for the products and feedstocks; eT=[1, 1, 1, 1, 1]; the mapping f in constraint (a) is a set of product quality models; qf are feedstock qualities; s are product specifications (see Tables 2 and 3); is a ‘‘bias’’ vector that is estimated from the available process measurements; constraint (b) represents the minimum and maximum limits for the product demands and feedstock availabilities; and A ¼ ½I55 e15 015 T B ¼ ½I55 015 e15 T bl ¼ 0:0 0:0 0:0 0:0 0:0
ð34Þ
583:33
666:67
T ð35Þ ð36Þ
bu ¼ ½1000:0 541:67 250:0 375:0 583:33 833:33 833:33 T
a
Deterministic solution is derived from the problem in which the random variables are replaced by their expectations.
ð37Þ
382
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
Fig. 2. Results for different numbers of Monte Carlo simulations. Table 2 Feedstock qualities Feedstocka
No. 1
No. 2
No. 3
No. 4
No. 5
RON MON Olefin Aromatics RVP
70.7 68.7 1.8 2.7 12.0
94.1 80.5 1.0 58.0 3.8
93.8 90.0 0.0 0.0 51.5
92.9 80.8 48.8 22.8 5.3
95.0 91.7 0.0 0.0 6.6
a Feedstock No. 1: LSR naththa; No. 2: reformate; No. 3: n-butane; No. 4; catalytic gas; No. 5: alkylate.
Fig. 3. Gasoline blending process.
Table 3 Product specifications Regular grade
The mixing rules used in the process model are the Ethyl RT-70 approach for RON and MON, and the Blending Index method for RVP (see Appendix A). Then, the process model f is a set of nonlinear equations, which is used to describe the ‘‘true’’ plant behaviour for the purposes of this simulation. In the RTO system, f is linearized about the nominal operating conditions (see Appendix B). The resulting plant/model mismatch is addressed using on-line bias update. The challenge of gasoline blending process is that there is some uncertainty in the qualities (e.g. RON, MON and
RON MON RVP
Premium grade
Min
Max
Min
Max
88.5 77.0 N/A
N/A N/A 10.8
91.5 80.0 N/A
N/A N/A 10.8
RVP) of each feedstock. In this case study, the quality disturbances are modeled as a Markov process [22]: q f ðtÞ ¼ q f ðt 1Þ þ C"ðtÞ
ð38Þ
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
Fig. 4. Reformate qualities.
Fig. 5. Blended product qualities—conventional RTO.
383
384
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
where: "ðtÞ follows the standard multivariate normal distribution; and C is chosen such that the standard deviation of the quality disturbance is 0.75% of nominal measurement value. The correlation matrix of the feedstock qualities (i.e. RON, MON and RVP) for each feedstock is: 2
1:0 R ¼ 4 0:85 0:0
3 0:85 0:0 1:0 0:0 5 0:0 1:0
ð39Þ
This correlation matrix would usually be determined via statistical analysis of laboratory data. Fig. 4 shows the variation in the reformate qualities used in the case study (the other feedstocks show similar, but not identical, variation in the case study). Note that the feedstock quality disturbances were assumed to remain unchanged during each RTO interval. Then, the optimization problem for this gasoline blending case study contains 10 decision variables (i.e. the flows of each of the 5 feedstocks to the regular and premium gasolines) and 15 uncertain parameters (i.e. the RON, MON and RVP of each of the 5 feedstocks).
5.1. Simulation results For the purpose of comparison, four RTO algorithms are used in this case study: (1) an ideal RTO, which provides the maximum possible profit from the blending process and meets all blended quality specification as well as limitations on product demands and feedstock availabilities, by using the perfect knowledge of the blending process (i.e. a perfect blending model) and all feedstock qualities, present and future. Although such an ideal RTO cannot be realized in practise, it can serve as a benchmark against which the performance of other RTO systems can be evaluated; (2) a conventional RTO, which includes a ‘‘bias’’ update for in the linearized process models using the available process measurements and then the solution of a linear program (LP) to determine the optimum operating conditions using the updated process model; (3) an IPC RTO T e x max cTp T cTf ðx þ yÞ x;y e y s:t: : P½fj ðx; y; qf Þ j 5 sj 5 j ðj ¼ 1; ; kÞ ðaÞ bl 4 Ax þ By 4 bu
ðbÞ ð40Þ
Fig. 6. Blended product qualities—IPC RTO.
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
Step 3.
and (4) a JPC RTO T e x cTf ðx þ yÞ max cTp T x;y e y s:t: : P½fðx; y; qf Þ 5 s 5 0
ðaÞ
bl 4 Ax þ By 4 bu
ðbÞ
ð41Þ
Step 4. In the simulations, the assumed probability levels in either IPC or JPC RTO are taken as 95% (i.e. 0 ¼ j ¼ 0:95; j ¼ 1; ; kÞ for all product quality constraints and the simulations were run for 100 RTO intervals. To evaluate the probability values and gradient vectors of the probability as a function of the decision variables, 5000 Monte Carlo simulations for each decision variable are carried out during every RTO interval. The simulation studies for the JPC RTO algorithm were carried out as follows: Step 1.
Step 2.
gather the process operating data (e.g. feedstock flow, blended gasoline qualities, etc.) once the process reaches a new steady state; linearize the process model at the current steady-state operating point and calculate the adjustable model parameters (bias terms) using the available process measurements;
385
solution of the JPC optimization problem: (a) solve problem (22) to find initial feasible solution; and solve the IPC optimization problem to find initial infeasible solution; (b) calculate the constraint probabilities and gradients using Monte Carlo simulations; (c) find the optimal solution using supporting hyperplane method; implement the new operation policy, wait for a new steady state and return to Step 1.
Fig. 5 shows the product quality for both gasoline grades determined using conventional RTO. As expected, the conventional RTO is not particularly robust to the uncertainty in the feedstock qualities. The conventional RTO produces frequent quality constraint violations in RVP for both gasoline grades and in MON for the premium grade. Further, these violations represent a substantial portion of the assumed feedstock quality variance. The simple IPC algorithm appears to considerably improve RTO robustness, even though it treats the probability constraints independently, as shown in Fig. 6. The frequency of constraint violation has been reduced by a factor of more than 3 (cf. Table 4); however, this constraint violation frequency is still substantial and may not be acceptable.
Fig. 7. Blended product qualities—JPC RTO.
386
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
Fig. 7 illustrates the performance of the JPC RTO and shows that only a few constraint violations are produced in the product qualities. Fig. 7 also shows that the constraint violation frequency for the JPC RTO is substantially less than either the conventional RTO or IPC RTO. Table 4 compares the conventional RTO, IPC RTO and JPC RTO approaches with respect to the number of quality constraint violations. It shows that number of constraint violations decreases from 76% for the Conventional RTO to 20–26% for the IPC algorithm, and to 6–12% by employing JPC algorithm. Note that the probability of constraint violations may exceed the specific probability levels for any given set of experiments due to the realization of the feedstock qualities used and plant/model mismatch.
Although the JPC algorithm improved the solution feasibility significantly, it did so at a considerable computational cost. The number of calculations required by the supporting hyperplane method and Monte Carlo simulations are considerable. Table 5 gives the time necessary to complete the entire simulation for each of the RTO approaches. Note that there was little difference between the computational costs of the conventional RTO and IPC RTO; whereas, JPC RTO required more than an order of magnitude increase in computation time beyond that required for the other methods. In the JPC case, the calculation time required at each RTO interval was approximately ten minutes. This was within the MATLAB1 environment and a substantial reduction in the required computing time could be expected in a production computing system. Table 5 RTO calculation times
Table 4 Number of product quality violations Method
Regular grade violations
Premium grade violations
Conventional RTO IPC RTO JPC RTO
76 26 6
77 20 12
Method
Time (s)a
Conventional RTO IPC RTO JPC RTO
20 25 6102
a All calculations were implemented in MATLAB1 version 5.3 (CPU: Pentium III/450 MHz; RAM: 64M; OS; Windows1 2000).
Fig. 8. Cumulative operating profit.
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
5.2. Economic performance analysis There are a number of difficulties in evaluating the operating profit for both conventional RTO and robust RTO, as the actual profit realized depends on how the offspecification product is valued. In general, off-specification product must be subjected to additional treatment. For this case study it was assumed that such additional treatment would increase the cost of production by 25% of the total product value. Fig. 8 gives the cumulative operating profit profiles for each RTO approach. The constraint violation penalty was arbitrarily chosen and affects the slope of the curves given in Fig. 8, but not the comparative ranking. Because of the off-specification product penalty, Conventional RTO gives a decreasing operating profit (or increasing operating loss). As expected, the JPC RTO outperforms all but the Ideal RTO. This is predominantly due to its low frequency of product quality violations. The IPC RTO approach produces slightly less operating profit than the JPC RTO approach, but proves to be a substantial improvement over the conventional RTO approach. Finally, it is worth noting that the distance between any two curves at a given RTO interval in Fig. 8 represents the difference in operating profit produced by the corresponding RTO algorithms. It can be concluded that the JPC approach garnered approximately 82% of the theoretically available operating profit within the simulation time span; whereas, the IPC RTO obtained approximately 64% of the theoretically available operating profit. This difference may provide the economic driving force for justifying the additional computing resources required to complete the JPC RTO calculations within an acceptable time frame for RTO purposes.
387
constrained problem. This increased computational load is primarily due to the Monte Carlo simulations required to characterize the joint probability density functions for the uncertain problem parameters and the computationally expensive optimization algorithm used for the JPC approach. Specially tailored algorithms, parallel processing computers, and probability distribution approximation methods may substantially reduce the computational requirements for such jointly constrained problems. The economic benefit derived from the additional robustness of the JPC RTO would then serve as a basis for economic justification of the required additional investment in computing resources. Alternatively, the individual constraint satisfaction probabilities in the IPC RTO could be ‘‘tuned’’ to provide a similar robustness to the JPC approach; however, such probability tuning is arbitrary and may be unmanageable for large optimization problems. Finally, a key issue that has not been properly addressed in the literature is that of dealing with nonlinear constraints in stochastic optimization. Since the process models used in the heart of many RTO systems are nonlinear, further work is required to develop methods to deal with such constraints within a robust RTO framework.
Acknowledgements The authors are grateful for the financial support of Imperial Oil Ltd., Sunoco Inc., and Natural Sciences and Engineering Research Council of Canada.
Appendix A. Gasoline blending model 6. Summary and conclusions Various sources of uncertainty, including market uncertainty, process uncertainty, measurement uncertainty and model uncertainty, exist in any RTO systems, and influence RTO performance. Then the challenge for any well-designed RTO system is to make robust decisions under such uncertainty in real-time. The work presented here developed a stochastic optimization approach to robust RTO that deals effectively with parametric uncertainty. The key contribution is the incorporation of uncertainty into the RTO problem using stochastic programming methods. Algorithms were presented that treated the probabilistic constraints either independently or jointly. When correlation exists among the probabilistic constraints, it was shown that ignoring this correlation reduced the robustness of the RTO system. Further, the increased computational cost required to solve the joint probability constrained problem was shown to be substantially more than that of the independently
RON model gRON ¼ ðeT xÞRON ðrT xÞðsT xÞ ¼ rT x þ 1 rT diagðsÞx ðeT xÞ ðoT xÞ2 þ 2 oTs x T ðe xÞ ðaT xÞ2 T þ 3 a s x T ðe xÞ
ð42Þ
where: x is the feedstock flow rate, r is the RON of each feedstock, RON is the blended RON, s ¼ r m, o is the olefin content of each feedstock, os is the Olefin content squared, a is the aromatics content of each feedstock, as is the square of the aromatics content, and i are the parameters in the model ða1 ¼ 0:03224; a2 ¼ 0:00101; a3 ¼ 0; a4 ¼ 0:0445; a5 ¼ 0:00081; a6 ¼ 0:00645Þ.
388
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
rrT gRON ¼ ðrT xÞxT þ ðsT xÞxT xT þ 1 xT diagð2r mÞ ðeT xÞ
MON model gMON ¼ ðeT xÞMON
ðmT xÞðsT xÞ ¼ mT x þ 4 mT diagðsÞx ðeT xÞ ðoT xÞ2 þ 5 oTs x T ðe xÞ 2 6 ðaT xÞ2 T a x T þ 10000ðeT xÞ s ðe xÞ
T rm gRON
rpT gRON ¼ 0T
gMON ¼ ðgMON Þ0 þ ðrxT gMON Þ0 ðx x0 Þ þ MON T þ ðrrT gMON Þ0 ðr r0 Þ þ ðrm gMON Þ0 ðm m0 Þ
þ ðrpT gMON Þ0 ðp p0 Þ
gRVPBI ¼ ðeT xÞRVPBI ¼ pT x
ð44Þ
where p is the RVP blending index for each feedstock and RVPBI is the blend index of RVP for the blended product.
rxT gMON ¼ mT 0 B T þ 4 B @m diagðsÞ 0 B T þ 5 B @os
RON linearized model
T þ ðrrT gRON Þ0 ðr r0 Þ þ ðrm gRON Þ0 ðm m0 Þ
ð45Þ
1 ðoT xÞðoT xÞ T e C ðeT xÞ C A T ðe xÞ
2ðoT xÞoT
12 ðaT xÞðaT xÞ T a x s 6 B ðeT xÞ C B C eT @ A T 10000 ðe xÞ 26 1 ðaT xÞðaT xÞ T a x ðeT xÞ 10000 ðeT xÞ s 0 1 ðaT xÞðaT xÞ T T T e 2ða xÞa B T C ðeT xÞ Ba C @ s A T ðe xÞ þ
where rxT gRON ¼ rT 0
1 ðrT xÞðsT xÞ ðr xÞs þ ðs xÞr BT ðeT xÞ C C þ 1 B r diagðsÞ @ A ðeT xÞ
B T þ 3 B @as
1 ðmT xÞðsT xÞ T e C ðeT xÞ C A T ðe xÞ
ðmT xÞsT þ ðsT xÞmT
0
gRON ¼ ðgRON Þ0 þ ðrxT gRON Þ0 ðx x0 Þ þ RON þ ðrpT gRON Þ0 ðp p0 Þ
ð50Þ
where
Appendix B. Model linearization
0
ð49Þ
MON linearized model
RVPBI model
B T þ 2 B @os
ð48Þ
ð43Þ
where m is the MON of each feedstock and MON is the blended MON.
0
T ðr xÞxT T x diagðrÞ ¼ 1 ðeT xÞ
ð47Þ
T
T
T
T
T
T
ðmT xÞxT rrT gMON ¼ 4 xT diagðmÞ ðeT xÞ
1
ðo xÞðo xÞ T e C ðeT xÞ C A ðeT xÞ
2ðoT xÞoT
1 ðaT xÞðaT xÞ T e C ðeT xÞ C A T ðe xÞ
2ðaT xÞaT
ð51Þ
ð52Þ
T rm gMON ¼ xT ðsT xÞxT ðmT xÞxT T þ 4 x diagðr 2mÞ ðeT xÞ
ð53Þ
rpT gMON ¼ 0T
ð54Þ
ð46Þ
Y. Zhang et al. / Journal of Process Control 12 (2002) 373–389
RVPBI linearized model gRVPBI ¼ðgRVPBI Þ0 þ ðrxT gRVPBI Þ0 ðx x0 Þ þ RVPBI þ ðrrT gRVPBI Þ0 ðr r0 Þ T gRVPBI Þ0 ðm m0 Þ þ ðrm
ð55Þ
þ ðrpT gRVPBI Þ0 ðp p0 Þ where rxT gRVPBI ¼ pT
ð56Þ
rrT gRVPBI ¼ 0T
ð57Þ
T gRVPBI ¼ 0T rm
ð58Þ
rPT gRVPBI ¼ xT
ð59Þ
References [1] D.C. White, Online optimization: what have we learned?, Hydrocarbon Processing 6 (1998) 55–59. [2] A. Kraslawski, Review of applications of various types of uncertainty in chemical engineering, Chem. Eng. Process. 26 (1989) 185–191. [3] M.G. Ierapetritou, J. Acevedo, E.N. Pistikopoulos, An optimization approach for process engineering problems under uncertainty, Computers Chem. Engng 20 (1996) 703–709. [4] R. Becker, S. Hall, B. Rustem, Robust optimal decisions with stochastic nonlinear economic systems, Journal of Economic Dynamics and Control 18 (1994) 125–147. [5] J.M. Mulvey, R.J. Vanderbei, Robust optimization of large-scale systems, Operations Research 43 (1995) 264–281. [6] J. Darlington, C.C. Pantelides, B.A. Tanyi, An algorithm for constrained nonlinear optimization under uncertainty, Automatica 35 (1999) 217–228.
389
[7] P. Kall, S.W. Wallace, Stochastic programming, Wiley, New York, 1994. [8] A. Pre´kopa, Stochastic programming, Kluwer Academic, Dordrecht, Netherlands, 1995. [9] J.F. Forbes, T.E. Marlin, Model accuracy for economic optimizing controllers: the bias update case, Ind. Eng. Chem. Res. 33 (1994) 1919–1929. [10] J.F. Forbes, T.E. Marlin, J.F. MacGregor, Model adequacy requirements for optimizing plant operations, Computers Chem. Engng. 18 (1994) 497–510. [11] Y. Zhang, D. Nadler, J.F. Forbes, Results analysis for trust constrained real-time optimization, Journal of Process Control 11 (2001) 329–341. [12] C. Loeblein, J.D. Perkins, B. Srinivasan, D. Bonvin, Economic performance analysis in the design of on-line batch optimization systems, Journal of Process Control 9 (1999) 61–78. [13] I. Miletic, T.E. Marlin, On-line statistical results analysis in realtime operations optimization, Ind. Eng. Chem. Res. 37 (1998) 3670–3684. [14] D.S. Monder, J.F. Forbes, Incorporation of Parametric Uncertainty in the Gasoline Blending Control Problem, presented at CSChE Meeting, 1998. [15] A.T. Schwarm, M. Nikolaou, Chance-constrained model predictive control, AIChE J. 45 (8) (1999) 1743–1752. [16] P. Li, M. Wendt, G. Wozny, Robust model predictive control under chance constraints, Computers Chem. Engng. 24 (2000) 829–834. [17] D.E. Kassmann, T.A. Badgwell, R.B. Hawkins, Robust steady state target calculation for model predictive control, AIChE J. 46 (2000) 1007–1024. [18] A. Charnes, W.W. Cooper, Deterministic equivalents for optimizing and satisfying under chance constraints, Operations Research 11 (1963) 18–39. [19] A. Pre´kopa, Numerical solution of probabilistic constrained programming problems, in: Y. Ermoliev, R.J.-B. Wets (Eds.), Numerical Techniques for Stochastic Optimization, Springer, Berlin, 1988, pp. 123–139. [20] T. Sza´ntai, A computer code for solution of probabilistic-constrained stochastic programming problems, in: Y. Ermoliev, R.J.B. Wets (Eds.), Numerical Techniques for Stochastic Optimization, Springer, Berlin, 1988, pp. 229–235. [21] F.Jr. Veinott, The supporting hyperplane method for unimodal programming, Operations Research 15 (1967) 147–152. [22] J. Honerkamp, Stochastic Dynamical Systems: Concepts, Numerical Methods, Data Analysis, VCH Publishers, New York, 1993.