Data-Driven Modeling and Optimization of Semibatch Reactors Using

In this study, a new data-driven approach has been proposed for modeling and trajectory optimization of a batch or a semibatch process. The approach i...
0 downloads 8 Views 224KB Size
Download PDF
Ind. Eng. Chem. Res. 2004, 43, 7539-7551

7539

Data-Driven Modeling and Optimization of Semibatch Reactors Using Artificial Neural Networks† K. Yamuna Rani* Process Dynamics and Control Group, Chemical Engineering Sciences, Indian Institute of Chemical Technology, Hyderabad 500 007, India

Sachin C. Patwardhan‡ Department of Chemical Engineering, Indian Institute of Technology, Bombay, Powai, Mumbai 400 076, India

In this study, a new data-driven approach has been proposed for modeling and trajectory optimization of a batch or a semibatch process. The approach is based on parametrization of input and output trajectories as finite-dimensional vectors using orthonormal polynomials (i.e., Fourier coefficients). Using input/output trajectory information available in historical databases, an artificial neural network (ANN) based model has been developed for capturing the dynamics of semibatch processes operated over a fixed interval of time. The parametrized input trajectories, initial states, and process parameters are considered as inputs to the ANN-based model, which predicts output trajectories in terms of Fourier coefficients. Single-rate as well as multirate systems can be modeled by this approach with equal ease. The resulting algebraic model is further used to formulate an optimal control problem, which can be solved using conventional nonlinear programming techniques to generate open-loop optimal input policies or optimal setpoint trajectories. The effectiveness of the proposed ANN-based modeling and trajectory optimization scheme is demonstrated using simulation studies on a benchmark multiple-input multiple-output semibatch process reported in the literature. Analysis of the simulation results reveals that the proposed ANN-based modeling approach is capable of capturing the nonlinear as well as the time-varying behavior inherent in the semibatch system fairly accurately. In addition, it also captures batch-to-batch variations in initial conditions and other process parameters. The results of operating trajectory optimization based on the proposed single-rate as well as the multirate ANN model are comparable to the results of trajectory optimization obtained using the exact first principles model. 1. Introduction There has been a tremendous growth in the chemical industry involving the use of batch and semibatch reactors because of the demand for a large number of specialty chemicals. Operation in a semibatch mode is typically used for strongly exothermic reactions because one can balance the reaction heat with the available cooling of the reactor and adjustment of the feed rate.1,2 Further, in biochemical reactions, fed-batch operation is preferred in cases where substrate overfeeding inhibits the growth of microorganisms. Because batch processing is predominantly used in the manufacture of low-volume and high-value products, even a marginal increase in the product yield can lead to a considerable improvement in profitability.3 Thus, model-based optimal control of batch and semibatch processes is becoming increasingly important for surviving under the competitive market conditions. The development of a mathematical model is the key step in any model-based optimization and control scheme. Models for physical systems can be classified as whitebox and data-driven models. White-box or mechanistic † IICT Communication Number 030611. * To whom correspondence should be addressed. Tel.: ++9140-27160123 ext 2287. Fax: ++91-40-27193626. E-mail: [email protected]. ‡ E-mail: [email protected].

models are constructed entirely from the first principles (laws of physics, thermodynamics, chemical kinetics, etc.), together with physical insight into the system behavior. The conventional approach to modeling of batch processes is based on white-box models. Most of the available optimization and model-based control schemes for batch/semibatch reactors have been formulated based on the mechanistic models.1 Recently, Srinivasan et al.4 have proposed the use of robust optimization in the absence of measurements and a measurements-based optimization approach when measurements are available in order to handle the problem of uncertainty. Even by this approach, a model based on first principles is necessary to formulate the nominal optimization problem. While these models are globally valid and are capable of capturing the system behavior quite accurately, the development of the first principles model is often a difficult and time-consuming exercise particularly when the system under consideration involves complex and multiple reactions. Moreover, all of the states necessary for the model development are typically not measurable. Therefore, from a practical viewpoint, simplified data-driven nonlinear models can provide an attractive alternative for formulating batch reactor optimization and control problems. Data-driven models can be classified as gray-box and black-box models depending on the extent of physical knowledge used in the model. Gray-box models are employed when some physical insight is available but

10.1021/ie0305521 CCC: $27.50 © 2004 American Chemical Society Published on Web 10/19/2004

7540

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004

some of the model components have black-box structure. During the past decade, artificial neural networks (ANNs) have gained importance as versatile data-driven tools for modeling nonlinear steady-state as well as dynamic processes.5 The most commonly employed ANN models for chemical processes have been the hybrid models, combining mass balance equations with ANN models for unknown kinetics.6 In the case of energy balance equations, it also becomes necessary to incorporate an estimate of heat released during the reaction into the hybrid models. Van Can et al.7,8 have developed a gray-box model combined with macroscopic balances for bioprocesses and have addressed issues on the types of data used, interpolation and extrapolation capabilities, etc. Costa et al.9 have employed ANNs combined with process theoretical knowledge for batch process optimization and control, whereas Schenker and Agrawal10 have used a combination of two ANNs with a mechanistic model in a model predictive control (MPC) framework. To develop (estimate model parameters) and employ gray-box models, it is required that the measurements of all of the state variables are available. In practice, however, some of the concentrations are not measurable online, and this leads to insufficient information in mass balance equations. Black-box models can be used when no physical insight is available, but the chosen model structure belongs to families that are known to have good flexibility and have been successful in the past.11-13 Studies on batch/semibatch process modeling suitable for optimization based on the process operating data without the knowledge of the underlying physical model are still at their infancy. Popular classes of black-box models are linear and nonlinear time-series models such as ARMAX and NARMAX models. Conventionally, ANNs have been used as nonlinear difference equations (NARX or NOE) models for capturing the dynamics of continuously operated processes.12 Zhang et al.14 proposed robust neural networks for predicting the polymer quality in batch polymerization reactors by stacking multiple nonperfect neural networks, which are developed based on the bootstrap resamples of the original training data. Recurrent neural networks have been employed in nonlinear autoregressive moving average (NARMA) models proposed for the prediction of the heat-transfer fluid temperature in batch reactors.15 Iatrou et al.16 proposed a scheme for modeling nonlinear nonstationary dynamic systems using ANNs with polynomial activation functions modulated by an appropriate time function. Chang and Hung17 have proposed the use of neural networks rate function models for optimization of batch polymerization reactors. The output derivatives with respect to time (generated by approximating the outputs as nonlinear functions of time) from a few typical batches are related to the output measurements and the temperature at every sampling instant with the help of a feed-forward neural network (FNN). The resulting model, combined with a proportional-integral controller, is used in the formulation of a modified two-step method for the determination of the optimal temperature trajectory. A peculiarity of batch and semibatch operations is the variations in the dynamics from batch to batch due to variations in the initial conditions, state of the processing equipment, etc. The development of ANNs as difference equation (or rate function) models, on the other hand, requires sequentially connected data, and

it is difficult to capture batch-to-batch variations through such time-series models. However, the ability of neural networks to learn a complex nonlinear relationship even when the input information is noisy and imprecise makes the neural network technology well suited for solving problems in batch process modeling. There have been some attempts to address the problem of batch-to-batch variations directly through batch-to-batch optimization methods18,19 based on datadriven models. The main aim in these approaches is to achieve the optimal solution iteratively over batches analogous to numerical optimization methods. Dong et al.18 have proposed an approach using neural networks for two semibatch reaction systems. The parametrized input variables are related to the state variables at the final time over batches with the help of a neural network multiway partial-least-squares algorithm. The optimization is performed by iterating from batch to batch with the help of the model obtained, where the gradients are evaluated analytically by incorporating information on the final state variable measurements of the previous batch. In this approach, the state variable trajectories are not a part of the model, and therefore it may not be possible to handle constraints on the state variables during the batch. Recently, Lee and Lee19 have proposed an approach combining iterative learning control and MPC for a semibatch reactor trajectory tracking and yield optimization. This approach is based on the representation of the process dynamics as a linear static map between the parametrized input and output trajectories in batch and regression models between the final product quality and the input/output trajectories. A batch-to-batch transition model is formulated with the help of an observer, and the input correction vector from batch to batch is computed using a MPC law. This approach, however, employs a linear representation of the input/output dynamics, and this restricts its applicability to batch processes with mildly nonlinear dynamics. In the present study, a novel modeling approach is proposed for the development of a black-box model for semibatch processes20,21 based on historical batch records, which combines some of the positive features of the above two approaches. The proposed modeling approach, based on parametrization of the input and output trajectories into finite-dimensional vectors, is discussed first. The ability of neural networks to learn complex nonlinear relationships has been exploited for the development of a black-box model. The problem of batchto-batch variations is handled indirectly in the present modeling approach by incorporating initial conditions and process parameters as additional inputs and is further consolidated by considering a large number of data sets for training and testing. This model is incorporated into a constrained optimization problem formulation, and the solutions are discussed. The efficacy of the proposed black-box modeling and optimization scheme is evaluated by conducting simulation studies on a multiple-input multiple-output (MIMO)22 semibatch process. 2. Data-Driven Modeling Approach In the present work, it is assumed that (i) the batch/ semibatch process under consideration can be adequately represented using lumped-parameter models, (ii) batch/semibatch process under consideration is controlled using a digital computer and the manipulated

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7541

inputs are piecewise constant, (iii) output trajectories and measured disturbances are available in the form of data vectors sampled at regular and sufficiently small time intervals, (iv) all of the batches are of identical time duration, and (v) a sufficiently large historical database is available for the development of data-driven models. Consider a general MIMO lumped-parameter process modeled by a vector differential equation represented as

dx ) f(x,u,d,p), x(t0) ) x0 dt

(1)

where x(t) ∈ Rn represents the system states, u(t) ∈ Rm represents the manipulated inputs, d(t) ∈ Rd represents the disturbances for t ∈ [t0, tf], p ∈ Rp represents the model parameters, and f represents an n-dimensional function vector. The system outputs, y ∈ Rr, are given by the relation

y(t) ) g[x(t)]

(2)

where g represents r-dimensional function vector mapping states to measured outputs. Let the state transition function χ[t;t0,x0,u(ut),d(ut),p] denote the solution of eq 1 obtained by implementing u(‚) with d(‚) over the set ut ≡ {τ: t0 e τ e t}such that

χ[t0;t0,x0,u(t0),d(t0),p] ) x0 x(t) ) χ[t;t0,x0,u(ut),d(ut),p] ) x0 +

∫ttf[x(τ),u(τ),d(τ),p] dτ 0

(3)

where the output at time t is given by eq 2. Let y(ut) represent the output profile over the time interval t0 e τ e t. Then,

y(ut) ) g{χ[ut;x0,u(ut),d(ut),p]} ) Ψ[y0,u(ut),d(ut),p]

(4)

where Ψ[‚] represents an algebraic map relating input and disturbance profiles to the output profile. In practice, it is seldom possible to develop such a closed-form map for a given system even if the first principles model is available. In the present work, it is proposed to develop an approximation to Ψ[‚] using input/output data. The concept of representing dynamic relationships of ordinary differential equations with initial conditions using “static” or “algebraic” maps is not entirely new to modeling and control literature. The classical transfer function representation is essentially an algebraic (static) map relating transformed inputs to transformed outputs. The transfer function maps infinite-dimensional vector u(s) to another infinite-dimensional vector y(s). The proposed modeling scheme can be viewed in an analogous manner, where the ANN represents an algebraic or static nonlinear map between the transformed input and transformed output signals. This analogy is schematically represented as shown in Figure 1. Similar to the transfer function map, the ANN map (as an approximation to Ψ[‚]) allows computation of the output profile for any arbitrary input profile u(t) and initial conditions. 2.1. Finite-Dimensional Approximation of Input and Output Trajectories. It should be noted that

Figure 1. Schematic representation of linear and ANN mapping analogy.

Ψ[‚] represents mapping from one set of continuous functions over t0 e t e tf (i.e., Cm[t0,tf] × Cd[t0,tf] × Rr × Rp) to another (i.e., Cr[t0,tf]). Here, C[t0,tf] represents the infinite-dimensional space of a real-valued continuous function over interval [t0, tf] with an infinite norm defined on it, i.e., for any f(t) ∈ C[t0,tf]

|f(t)|∞ ) max |f(t)| t∈[t0,tf]

(5)

The space of vectored continuous functions, such as Cm[t0,tf], can be constructed by “stacking” spaces C[t0,tf]. Because it is difficult to deal with infinite-dimensional vectors {u(t): t0 e t e tf}, {y(t): t0 e t e tf}, and {d(t): t0 e t e tf} directly in the development of an approximation for Ψ[‚], the first step in the model development is to approximate these infinite-dimensional vectors using low-order finite-dimensional vectors. We invoke the well-known Weierstrass approximation theorem for developing such a finite-dimensional representation. Weierstrass Approximation Theorem:23 The set of all polynomials with real coefficients is dense in the real space C[a,b]; i.e., given any  > 0, for every vector f(t) ∈ C[t0, tf], there exists a finite-order polynomial p(t) such that |f(t) - p(t)| <  for all t ∈ [a, b]. In other words, it is possible to represent every real continuous function as a finite-dimensional polynomial with reasonable accuracy. Therefore, in the present study, polynomial functions are chosen as the basis functions to generate an approximate finite-dimensional representation of input and output trajectories. Thus, the input trajectories, measured disturbance trajectories, and output trajectories are approximated as si

ui(t) ≈

aijφj(t) ∑ j)1

i ) 1, ..., m

ti

di(t) ≈

bijφj(t) ∑ j)1

i ) 1, ..., d

lk

yk(t) ≈

ckjφj(t) ∑ j)1

k ) 1, ..., r

(6)

where φj(t) represents the jth polynomial basis function, aij, bij, and ckj represent the polynomial coefficients for input, disturbance, and output trajectories, respectively, and si, ti, and lk denote the number of polynomial basis functions employed to approximate the input, disturbance, and output trajectories, respectively. This allows

7542

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004

representation of input and output trajectories in terms of finite-dimensional vectors, i.e.

Ui ≡ [ai1, ai2, ..., aisi]T = {ui(t): t0 e t e tf} Di ≡ [bi1, bi2, ..., biti]T = {di(t): t0 e t e tf} Yk ≡ [ck1, ck2, ..., cklk]T = {yk(t): t0 e t e tf}

φ1(t) ) e1(t); φ2(t) ) e2(t); φ3(t) ) e3(t); ... such that

(8)

where Y ) (Y1T, Y2T, ..., YrT)T, U ) (U1T, U2T, ..., UmT)T, and D ) (D1T, D2T, ..., DdT)T and Ω(‚) represents the ANN map that approximates the operator Ψ[‚]. Note that the independent inputs to the model include initial conditions and certain process parameters, which are generally constant throughout the batch duration. Unlike in continuous processes, the initial conditions play a significant role in batch/semibatch processes. Moreover, the optimal trajectories and performance also vary with the initial conditions. Therefore, it is necessary to incorporate the initial conditions into models used for optimization. Certain process parameters such as the catalyst concentration in a catalytic batch reaction or the reactant concentration in the feed stream in a semibatch reaction, etc., that remain approximately constant throughout the batch duration also influence the system optimal performance. Therefore, certain process parameters are also included as inputs to the models employed for optimization. The resulting ANNbased model can be trained using a historical database containing several batch or semibatch runs. 2.2. Approximation Using Orthonormal Polynomials. In practice, the most convenient norm to work with is the 2-norm

∫tt f(t)2 dt]1/2

|f(t)|2 ) [

f

0

(9)

because it is associated with an inner product. It can be easily shown that the Weierstrass approximation theorem also holds in L2[t0,tf], which represents completion of space C2[t0,tf] (i.e., the space of real-valued continuous functions with 2-norm defined on it as above).23 Thus, we make an additional assumption that {u(t): t0 e t e tf}, {y(t): t0 e t e tf}, and {d(t): t0 e t e tf} belong to the infinite-dimensional spaces Lm 2 [t0,tf], Lr2[t0,tf], and Ld2 [t0,tf], respectively. In other words, Ψ[‚] is now assumed to represent a map from Lm 2 [t0,tf] × Ld2 [t0,tf] × Rr × Rp to Lr2[t0,tf]. The main advantage of working with L2[t0,tf] is that it is a Hilbert space with an inner product defined on it as follows:

〈f(t), g(t)〉 )

∫tt f(t) g(t) dt f

0

(11)

(7)

Note that the above approach maps trajectories associated with entire batch or semibatch runs into a set of finite-dimensional vectors. The next step is to use these finite-dimensional representations in the development of an approximation for operator Ψ[‚]. We exploit the ability of ANNs (multilayer perceptron) to learn complex nonlinear algebraic relationships to develop such an approximation. Thus, an ANN map is constructed that relates input trajectories (U), disturbance trajectories (D), initial states (x0), and certain process parameters (p) to the output trajectory as

Y ) Ω(U,D,x0,p)

As a consequence, it is possible to define a set of orthogonal vectors in L2[t0,tf]. Also, the approximation problem can be posed as a least-squares estimation problem and can be solved analytically. Now, given an orthonormal polynomial basis for L2[t0,tf], say

(10)

〈ei(t), ej(t)〉 )

{

1 if i ) j 0 if i * j

(12)

any arbitrary function f(t) ∈ L2[t0,tf] can be expressed as a generalized Fourier series24

f(t) ) c1e1(t) + c2e2(t) + c3e3(t) + ... where

ci ) 〈f(t), ei(t)〉

(13)

and an n-dimensional best approximation in the leastsquares sense can be directly written as

f(x) = c1e1(t) + c2e2(t) + c3e3(t) + ... + cnen(t)

(14)

Significant advantages of using an orthonormal basis are (i) the approximation problem is well posed even for a high-order approximation and (ii) if the basis polynomials are chosen appropriately, the coefficients in the approximation decrease rapidly in magnitude with an increase in the polynomial order. Thus, only the first few coefficients are sufficient to generate a reasonably accurate approximation of a function. In the present work, we are dealing with continuous trajectories over the interval [t0, tf], which can be mapped to the interval [0, 1] by appropriate transformation of time coordinates. An orthonormal basis for L2[0,1] can be constructed by applying the GramSchmidt process25 to any linearly independent set of vectors belonging to L2[t0,tf], such as

φ1(t) ) 1; φ2(t) ) t; φ3(t) ) t2; φ3(t) ) t3; ... The resulting orthonormal polynomials are shifted Legendre polynomials (the first few are reported in Table 1). 2.2.1. Single-Rate Systems. There are two difficulties in applying the above approach to an approximation of input and output trajectories when the batch or semibatch process under consideration is controlled using a computer. In a computer-controlled process, the information about output trajectories (or measured disturbance trajectories) is available in the form of a sequence of values sampled at regular sampling intervals. The intersample behavior of the outputs (or measured disturbances) is not known. In addition, the sampled data are corrupted with measurement noise, and this calls for a different approach to the generation of the finite-dimensional approximation. On the other hand, the manipulated inputs are piecewise constant functions of time. In other words, the manipulated input trajectories are known over the entire time interval. Thus, we deal with input and output trajectories in a different manner while developing their finite-dimensional approximations. The true output trajectories are

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7543 Table 1. Orthonormal Polynomials over [0, 1] order, j

polynomial, ej

0 1 2 3 4 5 6 7 8 9 10

1 3.4641t - 1.7321 13.4164t2 - 13.4164t + 2.2361 52.9150t3 - 79.3725t2 + 31.7490t - 2.6458 210.0000t4 - 420.0000t3 + 270.0000t2 - 60.0000t + 3.0000 103(0.8358t5 - 2.0895t4 + 1.8573t3 - 0.6965t2 + 0.0995t - 0.0033) 104(0.3332t6 - 0.9995t5 + 1.1357t4 - 0.6057t3 + 0.1514t2 - 0.0151t + 0.0004) 104(1.3292t7 - 4.6522t6 + 6.4415t5 - 4.4733t4 + 1.6267t3 - 0.2928t2 + 0.0217t - 0.0004) 105(0.5306t8 - 2.1226t7 + 3.4668t6 - 2.9716t5 + 1.4286t4 - 0.3810t3 + 0.0519t2 - 0.0030t + 0.00004) 106(0.2119t9 - 0.9543t8 + 1.7976t7 - 1.8366t6 + 1.1030t5 - 0.3943t4 + 0.0810t3 - 0.0087t2 + 0.0004t - 0.00004) 106(0.2607t10 - 1.1020t9 + 1.8738t8 - 1.5940t7 + 0.6432t6 - 0.0332t5 - 0.0716t4 + 0.0268t3 - 0.0038t2 + 0.0002t - 0.00003) 106(0.0131t11 - 0.2897t10 + 1.0546t9 - 1.6304t8 + 1.2154t7 - 0.3271t6 - 0.1241t5 + 0.1188t4 - 0.0351t3 + 0.0046t2 - 0.0002t + 0.00003) 106(0.0030t12 + 0.0068t11 - 0.2885t10 + 1.0448t9 - 1.5718t8 + 1.1097t7 - 0.2304t6 - 0.1751t5 + 0.1348t4 0.0379t3 + 0.0049t2 - 0.0002t + 0.0003)

11 12

considered as smooth continuous functions, and measured data points are treated as sampled values of this function, which are corrupted with zero-mean white noise arising in the process of measurement. Thus, the coefficients of the polynomial approximation are determined as follows:

[ ][

e1(t0), ..., elk(t0) yk(t0) : : ) : : e1(tf), ..., elk(tf) yk(tf)

][ ]

ck1 : ) E‚c : cklk

(15)

The solution to this equation can be obtained in the least-squares sense as

c ) (ETE)-1‚ETy

(16)

The approximation for measured disturbances is developed in a similar manner. The input trajectories, on the other hand, are piecewise constant (over the sampling interval) and piecewise continuous functions of time over [t0, tf]. Because any piecewise constant function over the interval [t0, tf] can be arbitrarily closely represented by a continuous function, the piecewise constant (or continuous) functions belong to L2[t0,tf], i.e., completion of C2[t0,tf]. The input trajectories are approximated by a continuous function, which in the limit approaches a piecewise constant staircase-like function. A polynomial approximation of this continuous function can be obtained by taking the projections of the piecewise constant input trajectory in the direction of each orthonormal polynomial vector as follows: N

aij )

∑ (uik∫t k)1

tk+1

k

ej(τ) dτ) i ) 1, 2, ..., m; j ) 1, 2, ..., si (17)

where uik represents the ith input at the kth sampling instant and N represents the number of samples over the batch duration. Root-mean-square (rms) error and the Akaike information criterion (AIC) are employed to find the suitable order of polynomials for approximating input and output trajectories.26 The value of lk resulting in the minimum AIC with a reasonably small rms value is chosen as the order for the output trajectories. A similar procedure is

followed for the determination of order si for the input trajectories. 2.2.2. Multirate Systems. In the above development, it was assumed that the process outputs, manipulated inputs, and measured disturbances are sampled sufficiently fast and at identical sampling frequencies. In many batch or semibatch control applications, however, manipulated input moves are made at a much slower rate than the rate at which measurements are sampled.27 Also, in some applications, all of the measurements are not available at the same sampling frequency. Such batch or semibatch processes are, in fact, examples of multirate systems. In recent years, there has been growing interest in the dynamic modeling of multirate systems, and special methods are being developed for model identification using multirate sampled data.28-30 However, these developments are mostly limited to linear difference equation models developed in the neighborhood of a fixed operating point. The methods available in the literature for the development of datadriven NARX/NARMAX-type nonlinear difference equation models are not geared to deal with multirate sampled data systems. The proposed approach not only can deal with slowly sampled manipulated inputs but also can be easily adapted for systems where different outputs are sampled at different rates. When output measurements are available at different sampling rates, it is easy to see that eqs 15 and 16 can be used for each output without any modifications. In this section, we are concerned about modeling multirate systems where only a few manipulated input moves are made during the batch or semibatch operation while the measurements are collected at a sufficiently fast rate. It is assumed that all of the manipulated inputs are changed at identical frequencies, though this is not a strict requirement. For the multirate systems under consideration, the input trajectories are assumed to be piecewise constant over a time duration longer than the sampling interval, i.e.

ui(t) ) uki

for tk < t e tk + ∆ti

k ) 1, ..., si, i ) 1, ..., m

(18)

where ∆ti is defined as

∆ti ) tf/si

i ) 1, ..., m

(19)

7544

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004

In case tf is not an exact multiple of si, the length of all of the parts except the last one of the piecewise constant trajectory in terms of time is defined as

( )

f(x) )

tf ∆ti ) int si - 1

i ) 1, ..., m

(20)

∆tli ) tf - (si - 1)∆ti

i ) 1, ..., m

i ) 1, ..., m; j ) 1, ..., si

where

)0

)1

if

t + ∆ti t