Computers and Chemical Engineering 26 (2002) 863 – 873

Model decomposition based method for solving general dynamic

optimization problems

Yuanjie J. Huang, G.V. Reklaitis *, Venkat Venkatasubramanian

School of Chemical Engineering, Purdue Uni

6ersity, West Lafayette, IN

47907

-

1283

, USA

Received 24 July 2001; received in revised form 8 January 2002; accepted 16 January 2002

Abstract

The simultaneous strategy and control parameterization are the two most widely used direct methods for solving dynamic

optimization problems. The control parameterization approach generally results in a comparatively smaller nonlinear program
(NLP) but has difficulties in dealing with the path constraints. The advantage of the simultaneous approach lies in its ability to
handle path constraints, thus eliminating the need to obtain expensive and possibly infeasible intermediate solutions. The
disadvantage is that it requires solution of a potentially very large dimensional NLP. This work presents a decomposition strategy,
which combines the advantages of the control parameterization and simultaneous approaches for solving dynamic optimization
problems with path constraints. Using the proposed strategy, first the set of state variables x is divided into two sets x

1

and x

2

,

and the system is partitioned into two corresponding sub-systems. The criterion used to partition the state variables and the
system model are that the equations which define the state variables involved in the path constraints should be in the same
sub-system. For the resulting sub-systems, one is enforced in the master NLP through collocation method as in simultaneous
approach, the other is solved together with the sensitivities in a differential and algebraic equations (DAE) solver. This strategy
constitutes a general approach in the sense that for problems with a specific structure the method is equivalent to the control
parameterization method, while for other problems with special structures the approach is the same as the simultaneous approach.
In general it possesses the advantage of the simultaneous method in handling the path constraints since it is directly enforced in
the master NLP as well as the advantage of control parameterization in resulting a small master NLP because only a fraction of
the state variables are directly discretized. The method is demonstrated through several numerical examples. © 2002 Elsevier
Science Ltd. All rights reserved.

Keywords

:

Nonlinear programming; Dynamic optimization; Differential and algebraic equations; Collocation method

www.elsevier.com/locate/compchemeng

1. Introduction

Dynamic optimization problems often arise in pro-

cess engineering. Typical examples include finding an
optimal trajectory to perform the state transition task,
or finding a temperature profile that maximize the
selectivity to the desired product in a batch reactor.
Formulation of such problems results in a dynamic
optimization problem with path constraints on state
variables. The dynamic process model is composed of

sets of differential and algebraic equations (DAE). The
DAE system consists of differential equations that de-
scribe the behavior of the system and algebraic equa-
tions that represent the thermodynamic and kinetic
relationships of the system. Path constraints include
hard equipment constraints and operational con-
straints. Hard constraints represent physical limitations
on equipment, e.g. valves fully opened or closed. Oper-
ational constraints represent operating limits arising
from economic, safety or environmental considerations.

The solution algorithms for dynamic optimization

problems can be classified into direct methods and
indirect methods. Indirect methods formulate the prob-
lem in variational form and express the necessary con-
ditions for optimality as a two-point boundary value
problem (TPBVP), then use various iterative methods
to obtain the solution of the resulting TPBVP, as

* Corresponding author. Tel.: + 1-765-494-4075; fax: + 1-765-494-

0805.

E-mail

addresses

:

[e-mail]

(Y.J.

Huang),

[e-mail]

(G.V. Reklaitis),

[e-mail]

(V. Venkatasubramanian).

0098-1354/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved.
PII: S 0 0 9 8 - 1 3 5 4 ( 0 2 ) 0 0 0 0 7 - 8

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

864

shown in Denham and Bryson (1964), Dreyfus (1962),
Dixon and Bartholomew Biggs (1981), Miele, Mohanty,
Venkataraman and Kuo (1982). Using direct methods
the infinite dimensional dynamic optimization problem
is transformed into a finite dimensional nonlinear pro-
gram (NLP) and then solved using any NLP solver.
Direct methods are easier to implement since there is no
need to generate the additional costate/adjoint equa-
tions. In addition TPBVPs, especially those of large
dimensionality, are extremely difficult to solve numeri-
cally even with current computers.

The most widely used direct methods include control

parameterization techniques and simultaneous strategy
approaches. In the control parameterization techniques,
the control variables are discretized over finite elements
using polynomials. The coefficients of the polynomial
and the mesh size are the decision variables in a master
NLP (Brusch & Schapelle, 1973; Ray, 1981; Morison,
1984; Vassiliadis, 1993). State profiles and sensitivities
are evaluated through a separate DAE solver. In a
more recent work (Feehery, 1998), the author pointed
out that equality path constrained dynamic optimiza-
tion problems are essentially high-index DAEs and
developed dummy derivative method to solve these
high-index DAEs. The author also shown that inequal-
ity constrained dynamic optimization problems are es-
sentially

hybrid

discrete/continuous

systems.

This

hybrid characteristic arises because of the need to make
decisions concerning the ordering of constraint events
during the solution process. This combinatorial issue in
the solution of inequality constrained dynamic opti-
mization problems is one of the reasons that make such
problem difficult to solve. By contrast, in the simulta-
neous strategy method both the control and the state
variables are discretized using monomial representation
for piecewise polynomials on finite elements (Ascher,
Mattheij & Russell, 1988), and the coefficients and
interval lengths become the decision variables in a
larger NLP (Bartholomew-Biggs, 1995; Neuman & Sen,
1973; Vasantharajan & Biegler, 1990; Tanartkit &
Biegler, 1995).

Direct methods have been demonstrated to work

reliably for very large problems. Both methods have
their own advantages and disadvantages. Compara-
tively, the control parameterization method is easier to
use on large problems, since the system model equa-
tions are enforced through a numerical IVP integrator,
rather than including them as a set of constraints in the
master NLP. Thus, the resulting NLP is fairly small in
dimension. The main difficulty of this method is in
dealing with the path constraints: the equality path
constraints will cause the problem to be a high-index
DAE and the inequality path constraints will introduce
combinatorial issues in the solution process. The advan-
tages of the simultaneous approach lie in its treatment
of the path constraints as well as eliminating the need

to obtain expensive and possibly infeasible intermediate
solutions. The main shortcoming of this method may
be the large dimensionality of the resulting NLP.

In this paper, we present a model decomposition

strategy for solving dynamic optimization problem with
path constraints on state variables. In this method, the
system dynamic model is decomposed into two subsys-
tems, one is enforced in the master NLP together with
the path constraints through collocation methods, and
the other is satisfied through direct solution in a numer-
ical IVP integrator. The strategy seeks to combine the
control parameterization approach and simultaneous
approach by using the simultaneous strategy for the
first subsystem and solving the second subsystem in a
DAE solver. With control parameterization and simul-
taneous approach as its two special cases, the decompo-
sition strategy is a general direct method in this sense.
Therefore, it possesses both the advantages of simulta-
neous method and control parameterization method. In
addition, although a general method which can be used
for any dynamic optimization problems, it is extremely
well suited for some types of large scale problems that
often arise in process engineering.

In the next section, the details of the decomposition

strategy and some special cases are discussed. Then the
discretization strategies and sensitivity evaluation al-
gorithms and the DAE solver used in this work are
briefly described. Following this, several examples are
solved using the proposed decomposition method. The
paper ends with some conclusions and remarks.

2. Decomposition strategy

Consider the following general dynamic optimization

problem with path constraints:

(P1)

min

u(t),t

f

J =

„(x(t

f

), t

f

) +

&

t

f

t

0

L(x, u, t)dt

(1)

subject to the DAE:

f(x

; , x, u, t)=0

(2)

g(x

; , x, u, t)50

(3)

ƒ(x;(t

0

), x(t

0

), t

0

) = 0

(4)

where x

R

n

x

is state vector and u

R

n

u

is control vec-

tor. The dynamic process model (Eq. (2)) is assumed to
be a DAE with index

51. High-index DAE and equal-

ity path constrained system can be transformed into an
equivalent DAE with index equal to 1 or less using the
dummy derivative method as described in Feehery
(1998), Mattsson and So¨derlind (1993). Inequality (Eq.
(3)) is the path constraint.

To solve problem (P1), the large dimension of the

NLP resulting from a simultaneous method and the
difficulty in dealing with the path constraints when

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

865

using a control parameterization method motivated the
development of a new method which will overcome
these drawbacks. The goal of the decomposition strat-
egy is to ease the handling of path constraints while
minimizing the number of variables that need to be
discretized to reduce the size of the resulting NLP. For
the given problem (P1), all state variables that appear
in the path constraint (Eq. (3)) are relabeled x

1

R

n

1

.

The rest of the state variables are labeled as x

2

R

n

2

,

where n

x

= n

1

+ n

2

. Thus, the path constraint is in the

following form:

g(x

;

1

, x

1

, u, t)

50

(5)

With this partition of the state variables, the system

model (Eq. (2)) and the initial condition (Eq. (4)) are
then accordingly partitioned into:

f

1

(x

;

1

, x

;

2

, x

1

, x

2

, u, t) = 0

(6)

ƒ

1

(x

;

1

(t

0

), x

;

2

(t

0

), x

1

(t

0

), x

2

(t

0

), t

0

) = 0

(7)

and

f

2

(x

;

1

, x

;

2

, x

1

, x

2

, u, t) = 0

(8)

ƒ

2

(x

;

1

(t

0

), x

;

2

(t

0

), x

1

(t

0

), x

2

(t

0

), t

0

) = 0

(9)

such that matrices

#f

1

/

#q

1

,

#f

2

/

#q

2

,

#ƒ

1

/

#r

1

and

#ƒ

2

/

#r

2

are all nonsingular, where q

1

, r

1

, q

2

and r

2

are the

highest order time derivatives of x

1

and x

2

in Eqs.

(6) – (9), respectively. This is a sufficient condition for
Eqs. (6) and (8) to have an index (Brenan, Campbell &
Petzold, 1996) equal to or less than unity.

The general idea of the decomposition algorithm is as

follows. The first part of the model Eqs. (6) and (7) and
the path constraints (Eq. (5)) become the model equa-
tions and path constraints in a new dynamic optimiza-

tion

problem

which

is

solved

by

employing

a

simultaneous strategy. Both the input vector u and the
state vector x

1

are discretized and an NLP problem is

obtained. The profiles of x

2

and the sensitivities needed

by the NLP solver are directly obtained from a separate
DAE integrator, similar to a control parameterization
approach. The diagram of the decomposition algorithm
for solving general dynamic optimization problems is
illustrated in Fig. 1.

To partition the model equation f into two sub-sys-

tems f

1

and f

2

, one way is to obtain a maximum

transversal of the equation f and variable x. The maxi-
mum transversal gives a maximum or perfect matching
between the equations and the variables. Thus, function
f can be partitioned accordingly for known x

1

and x

2

.

An efficient algorithm, MC21A, for finding an aug-
menting path or maximum transversal for a given set of
equations and variables was described in Duff (1981).
Algorithm MC21A is adequate for our purpose.

Since functions f, f

1

and f

2

are all index-1 DAEs, by

the implicit function theory, there exists function f

0

2

such that x

2

= f

0

2

(u, t). Thus the subsystem (Eq. (6)) and

the objective function (Eq. (1)) can be written as
follows:

f

0

1

(x

;

1

, x

1

, u, t) = 0

(10)

J =

„

1

(x

1

(t

f

), t

f

) +

&

t

f

t

0

L

1

(x

1

, u, t)dt

(11)

With these modifications the original problem (P1) is

transformed into a dynamic optimization problem (P2)
as shown in Eqs. (12) – (15). Generally, problem (P2) is
much smaller in size than the original formulation (P1):

(P2)

min

u(t),t

f

J =

„

1

(x

1

(t

f

), t

f

) +

&

t

f

t

0

L

1

(x

1

, u, t)dt

(12)

subject to the DAE:

f

0

1

(x

;

1

, x

1

, u, t) = 0

(13)

g(x

;

1

, x

1

, u, t)

50

(14)

ƒ

1

(x

;

1

(t

0

), x

1

(t

0

), t

0

) = 0

(15)

Problem (P2) is then transformed into the following

NLP by employing a simultaneous approach (Vasan-
tharajan & Biegler, 1990; Logsdon, 1990; Cervantes &
Biegler, 1998):

min

z

P(z)

(16)

subject to:

h(z) = 0

(17)

z

L

5z5z

U

(18)

where z = (z

x

, z

u

, z

p

)

T

, z

x

and z

u

result from the dis-

cretization of x

1

and u in Eqs. (12) – (15), respectively,

and z

p

is the collection of all other variables. The NLP

solver needs the value of the gradients:

#h/#z

x

,

#h/#z

u

Fig. 1. The decomposition algorithm.

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

866

and

#h/#z

p

. Since Eq. (17) is derived from Eqs. (13) and

(14), derivative

#h/#z

x

does not depend on the subsys-

tem (Eqs. (8) and (9)). Thus z

x

does not need to be

included in the parameter set with respect to which the
sensitivity equations of Eq. (8) are obtained and solved.
This guarantees that the number of total equations
integrated in the DAE solver will not exceed the num-
ber of equations integrated when a control parameteri-
zation method is directly used to solve problem (P1).

The state variable x

1

takes different roles in Eqs. (2)

and (8) since x

1

in Eq. (2) is a dependent state variable

while x

1

in Eq. (8) can be considered as an input

variable to the sub-system consists of Eqs. (8) and (9).
Accordingly, we have two different methods to solve
the original dynamic optimization problems with re-
spect to whether x

1

is directly solved in the DAE

integrator or not. When x

1

is considered as an input to

Eq. (8), then only Eq. (8) and its corresponding sensi-
tivity equations need to be integrated in an IVP solver.
In the case when x

1

is deemed as a dependent state

variable, in order to obtain the profile of x

2

and the

sensitivities

#x

2

/

#p, the original model and its sensitivity

equations should also be integrated. In either case, the
number of equations (excluding the sensitivity equa-
tions) to be integrated is less than the number of
equations to be integrated in the control parameteriza-
tion method.

In solving the resulting NLP a good initial guess for

the decision variables is necessary. Rather than specify
the initial values for all the variables a better way to do
this is to specify the initial values for the control
variable u only and then to solve the original DAE (Eq.
(2)) to obtain the initial value for x

1

.

When the model (Eqs. (2) – (4)) in problem (P1) is in

special form, the decomposition method may reduce to
the

well-known

simultaneous

method

or

control

parameterization method as in the following two ex-
treme cases.

“

Control parameterization: only control variables are
discretized in this method, which is the case for the
decomposition method when there is no state vari-
ables in path constraint (Eq. (3)), the constraints
contain only control variables.

“

Simultaneous strategy: both state and control vari-
ables are discretized in the method, which is the case
for the decomposition approach when all state vari-
ables appear in constraints (Eq. (3)).

Let n

d

be the number of state variables which are

discretized in the method, then in the simultaneous
strategy, n

d

= n

x

; in the control parameterization

method, n

d

= 0; in the decomposition method, n

d

lies

somewhere between 0 and n

x

, and n

d

= n

x

or 0 are the

two special cases discussed above. Compared with the
simultaneous approach, the decomposition method dis-
cretizes a smaller number of state variable, thus result-
ing in a smaller size NLP than in the simultaneous

approach. Since all path constraints are included in the
first sub-system that is solved by employing a simulta-
neous method, path constraints handling becomes
much easier than in the control parameterization
method. The intermediate solution from the DAE
solver for a given input u is always feasible since u is a
feasible input. Thus, one can avoid the potential expen-
sive integrations only to find that the solution is infeasi-
ble. Most of the path constraints contained in the
dynamic optimization problems arising in process engi-
neering are in the form of lower or upper bounds on
the state variables which are directly transformed into
lower or upper bounds on the decision variables in the
resulting NLP. Such nonlinear constraints become lin-
ear constraints.

3. Discretization strategies

The decomposition method approximates the control

profiles using polynomials over a set of N finite ele-
ments of length h

i

= t

i + 1

− t

i

, i = 1, …, N. For the

convenience of bounding the control trajectories, the
controls are parameterized using either piecewise con-
stant functions when the control profiles are discontinu-
ous or Lagrange polynomials when the profiles are
continuous. For control variable u

j

in interval i the

Lagrange approximation is:

u

j

i

(

q)= %

K

k = 1

u

ijk

L

k

(K)

(

q) q[0, 1], i=1, …, n

u

(19)

where:

L

k

(K)

(

q)= 5

K + 1

l = 1,l

"k

q−q

l

q

k

−

q

l

(20)

and

q is the normalized time in the interval i:

q=

t − t

i

t

i + 1

− t

i

(21)

Those state variables that need to be discretized are

approximated as follows (Huang, Reklaitis & Venkata-
subramanian, 2000):

x(t) =

%

3

k = 0

c

k

i



t − t

i

h

i



k

,

for t

[t

i

, t

i + 1

], i = 1, …, N

(22)

with

c

0

i

= x

i

(23)

c

1

i

= h

i

f

i

(24)

c

2

i

= − 3x

i

− 2h

i

f

i

+ 3x

i + 1

− h

i

f

i + 1

(25)

c

3

i

= 2x

i

+ h

i

f

i

− 2x

i + 1

+ h

i

f

i + 1

(26)

where f

i

= x

; (t

i

), h

i

= t

i + 1

− t

i

. In order to derive Eqs.

(23) – (26) the state profile should be continuous and its

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

867

first order derivative also continuous. The continuity of
the state variables is reasonable because the differential
variables represent physical conserved quantities like
mass, energy or momentum. But for the derivatives this
condition is quite restrictive. A more practical condi-
tion is that on each interval [t

i

, t

i + 1

] the first order

derivative of x is only right continuous at t

i

and left

continuous at t

i + 1

. Under this condition, Eqs. (23) –

(26) can still be used when

f

i

= x

; (t

i

)

(27)

and

f

i + 1

= x

; (t

i + 1

)

(28)

on the ith interval, i = 1, …, N are replaced by

f

i

= x

; (t

i

+

)

b lim

t

“t

i

+

x

; (t)

(29)

f

i + 1

= x

; (t

i + 1

−

)

b lim

t

“t

i + 1

−

x

; (t)

(30)

4. Sensitivity evaluation

Sensitivity information must be provided for any

gradient-based NLP solver. Since both the objective
function and constraints can be formulated as model
equations, only the sensitivities of Eq. (2) are consid-
ered here. Taking derivatives with respect to the opti-
mization

parameters

p

results

in

the

following

sensitivity equation:

#f

(x;

s

; +

(f

(x

s +

(f

(u

(u
(p

= 0

(31)

with the sensitivity variables:

s =

(x
(p

(32)

The methods used to solve the combined DAE (Eq.

(2)) and its sensitivity system (Eq. (31)) can be classified
into two categories.

“

The staggered direct scheme in which the linear
systems for the sensitivity corrector steps are solved
directly after the convergence of the nonlinear cor-
rector step (Caracotsios & Stewart, 1985; Haug &
Ehle, 1982; Rabitz, Kramer & Dacol, 1983).

“

The simultaneous corrector method solves the DAE
and its sensitivities simultaneously, obtaining the
solution from one corrector iteration on the com-
bined nonlinear system (Maly & Petzold, 1996) or
two correctors iteration, one for the DAE and the
other for the sensitivities (Feehery, Tolsma & Bar-
ton, 1997).

From the computational point of view, the simulta-

neous methods are better than the staggered schemes
because the number of matrix factorizations that need
to be performed along the solution trajectory is mini-
mized. In this work, the DAE and sensitivity solver

DASPKSO (Maly & Petzold, 1996), an extension of the
DAE solver DASPK (Brenan et al., 1996) through the
implementation of the simultaneous corrector method
and use of an adaptive difference directional derivative
approximation to the sensitivity equations, is used to
calculate the state profiles and the sensitivities.

In the case that the input u is discontinuous at some

t

j

, 1

BjBN, then at t

j

a new state value instead of the

current state should be used as the initial value for the
DAE solver to continue integration in the time horizon
since some algebraic state variables may also have
discontinuities at t

j

.

5. Numerical examples

5

.

1

. Catalyst mixing problem

This problem determines the optimal mixing policy

of two catalysts along the length of a tubular reactor.
The mixing ratio of the catalysts is the control variable
and its optimal profile has a singular segment, making
this a difficult problem to solve. This problem was first
introduced in Gunn and Thomas (1965) and has been
solved by various researchers (Vassiliadis, 1993; Logs-
don, 1990; Tanartkit & Biegler, 1997; Bell & Sargent,
2000). The dynamic optimization formulation is de-
scribed as follows:

max

u

x

3

(1.0)

(33)

subject to:

x

;

1

= u(10x

2

− x

1

)

(34)

x

;

2

= u(x

1

− 10x

2

) − (1 − u)x

2

(35)

x

3

= 1 − x

1

− x

2

(36)

x(0) = [1.0

0

0]

(37)

u(t)

[0 1.0]

(38)

This example is used here because the dynamic opti-

mization formulation of this problem does not contain
any path constraint for state variables. Thus from our
model decomposition based strategy point of view, this
case is an extreme. The numerical algorithm of solving
this problem reduces to the control parameterization
method. Since there are no state variables involved in
path constraint in the formulation, the discretization of
the control variable only leads to an NLP with upper
and lower bounds on decision variables and no other
type of constraints on the decision variables. The di-
mension of this unconstrained NLP is proportional to
the number of mesh intervals used in the discretization.

Discussion on the details of applying control parame-

terization method to this problem is not the purpose of
this paper. For those interested, these can be found in
Vassiliadis (1993).

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

868

Table 1
Computational summary for fed-batch penicillin fermentation

NLP tolerance

DAE tolerance

Number of

Objective
function

iterations

l.E−7

−1.1228E6

32

l.E−7

process model into two corresponding sub-systems. The
incident matrix of Eqs. (40) – (44) is

X

P

V

v q

40

1

1

1

41

1

1

1

1

42

1

1

43

1

44

1

(48)

Using the algorithm MC21A, the following matching

is obtained:

X

P

V

v

q

40

×

1

1

41

1

×

1

1

42

1

×

43

×

44

×

(49)

where the × elements are the maximum transversal.
Thus, Eq. (40) that solves for variable X and the path
constraint, inequality (Eq. (46)) are enforced in the
master NLP. The rest of the system model equations
together with the sensitivity equations which constitute
a DAE system are solved using DASPKSO. In this
case, n

u

= n

1

= 1, n

x

= 5. Thus, if the same discretization

scheme is used, the number of variables of the NLP
resulting from decomposition method is about (n

u

+

n

1

)/(n

u

+ n

x

)

:30% of these arising in the simultaneous

method. The solution statistics are shown in Table 1.
The optimal control profile and some selected state
variable trajectories are shown in Figs. 2 and 3. These
results are very close to those reported in Feehery et al.
(1997); San and Stephanopoulos (1989).

5

.

3

. Pressure-constrained batch reactor

The reactions taking place in the reactor are:

A

X

k

1

k

2

2B

(50)

A + B

“

k

3

D

(51)

The dynamic optimization problem is described as

follows:

min

F

J = C

D

(t

f

)

(52)

subject to:

C

:

A

= − k

1

C

A

+ k

2

C

B

C

B

+

F

V

− k

3

C

A

C

B

(53)

C

:

B

= k

1

C

A

− k

2

C

B

C

B

− k

3

C

A

C

B

(54)

C

:

D

= k

3

C

A

C

B

(55)

5

.

2

. Fed-batch penicillin fermentation

The determination of the optimal feeding profile of

fed-batch fermentation requires the solution of a singu-
lar optimal control problem. The problem is considered
difficult due to the nonlinear dynamics of the system
model, the singular nature of the problem and the
existence of constraints on both the state and control
variables. The model used here was described and
solved in San and Stephanopoulos (1989). Following is
the problem formulation:

min

S(t)

J =

&

t

f

t

0

( −

qXV+0.0103PV+0.0744vX

+ 0.00102XV + 6913.58)dt

(39)

subject to:

X

: =vX−

FX

V

(40)

P

: =qX−0.01P−

FP

V

(41)

v=

0.11S

S + 0.006X

(42)

q=

0.004

1 + 0.0001/S + S/0.1

(43)

V

: =F

(44)

0.001

5S50.5

(45)

X

541

(46)

[X(0), P(0), V(0)] = [1.0, 0, 2.5E5]

(47)

The performance index in Eq. (39) represents the

overall profit. The biomass concentration X (g/l), the
amount of existing penicillin product P (activity per l),
the specific growth rate

v, the specific product forma-

tion rate

q and the reactor volume V (l) are the state

variables in this system. The substrate concentration S
(g/l) is the control variable. The feed rate F was set to
1666.67 l/h.

In this formulation, there is only one path constraint

that only contains one state variable X. Thus according
to the decomposition strategy described in Section 2,
among three state variables, only X is discretized over
finite elements using polynomials. Also control variable
S is approximated using piecewise constants. Then
Duff’s algorithm, MC21A, is used to partition the

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

869

Fig. 2. The optimal control profile.

Fig. 3. The optimal state trajectories.

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

870

Table 2
Computational summary for batch reactor

NLP tolerance

Objective

DAE tolerance

Number of
iterations

function

l.E−6

11.7446

562

l.E−6

C

A

C

B

C

D

N

P

53

1

1

54

1

1

55

1

1

1

56

1

1

1

1

57

1

1

(61)

the following match is obtained

C

A

C

B

C

D

N

P

53

×

1

54

1

×

55

1

1

×

56

1

1

1

×

57

1

×

(62)

The result of MC21A shows that Eq. (57) is to be

used to determine the constrained state variable P.
Thus, according to the decomposition strategy, Eq. (57)
and the path constraint Eq. (58) are enforced in the
master NLP. The rest of the system model and sensitiv-
ity equations are solved in the DAE solver. In this case,
n

u

= n

1

= 1, n

x

= 5. The NLP solved is again about 30%

of the size of the NLP solved in a simultaneous
method. In solving this problem the time interval, [t

0

,

t

f

], is discretized into 11 intervals. There are 35 vari-

N = V(C

A

+ C

B

+ C

D

)

(56)

PV = NRT

(57)

P

5340 000

(58)

0

5F58.5

(59)

[C

A

(0), C

B

(0), C

D

(0)] = [100, 0, 0]

(60)

This problem was also solved in Feehery (1998). The

rate constants are k

1

= 0.8 per h, k

2

= 0.02 m

3

/(mol h),

k

3

= 0.003 m

3

/(mol h), the volume V = 1.0 m

3

, the

temperature T = 400 K. As in the previous example,
there is only one path constraint on the state variable,
P. Thus P is the only state variable that should be
discretized. Applying Duff’s algorithm MC21A to this
model’s incident matrix:

Fig. 4. The optimal flowrate profile.

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

871

Fig. 5. The optimal pressure trajectory.

Fig. 6. The optimal concentration trajectories.

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

872

ables

and

23

equality

constraints

in

the

NLP

formulation. The computation statistics are given in
Table 2. Figs. 4 – 6 show the optimal profiles of the
control and some selected state variables. As shown in
the solution profiles, during the time horizon that the
state constraint is active the optimal control is highly
nonlinear which makes this problem difficult for
general control parameterization methods.

6. Conclusion

This paper has presented a scheme for solving dy-

namic optimization problems with path constraints on
state variables. The method is based on a model decom-
position strategy. For given path constraints, first the
state variables x are divided into two sets x

1

and x

2

. All

control variables and x

1

are then discretized using any

discretization strategy, such as the one described in
Section 3. All the path constraints and those model
equations that define vector x

1

are enforced in the

master NLP as constraints, while objective and the
sensitivities with respect to the control variables are
solved in a separate DAE and sensitivity solver. The
index of the DAE remains the same from iteration-to-
iteration since the profiles obtained from the integrator
are always feasible. In the case in which one of the
variable sets x

1

or x

2

is empty, this approach reduces to

the well known simultaneous method or the control
parameterization method. In this sense the decomposi-
tion strategy is a general method.

Using the proposed decomposition method, the di-

mension of the resulting NLP problem is in general
smaller than the one arising in the simultaneous ap-
proach. In most real cases, n

u

n

x

, and n

1

is about

10 – 30% of n

x

. If the proposed discretization strategy is

used, then the number of constraints in the NLP result-
ing arising in the decomposition approach is only (n

u

+

n

i

)/(n

u

+ n

x

), or 10 – 30% of the number from the

simultaneous approach. In terms of path constraints
handling, the decomposition method is preferable to the
control parameterization approach since a simultaneous
type of approach can be used, which can accommodate
the path constraints easily. Most of the path constraints
on the state variables in real cases appear as lower or
upper bounds, which become the lower or upper
bounds

on

the

decision

variables

in

the

NLP

formulation.

References

Ascher, U. M., Mattheij, R. M., & Russell, R. D. (1988). Numerical

solution of boundary

6alue problems for ordinary differential equa-

tions. Englewood Cliffs, NJ: Prentice Hall.

Bartholomew-Biggs, M. C. (1995). A penalty method for point and

path state constraints in trajectory optimization. Optimal Control
Applications and Methods,

16

(4), 291 – 297.

Bell, M. L., & Sargent, R. W. H. (2000). Optimal control of inequal-

ity constrained dae systems. Computers and Chemical Engineering,

24

(11), 2385 – 2404.

Brenan, K. E., Campbell, S. L., & Petzold, L. R. (1996). Numerical

solution of initial-

6alue problems in differential–algebraic equa-

tions. Philadelphia, PA: SIAM.

Brusch, R., & Schapelle, R. (1973). Solution of highly constrained

optimal control problems using nonlinear programming. AIAA
Journal,

11

, 135 – 136.

Caracotsios, M., & Stewart, W. E. (1985). Sensitivity analysis of

initial value problems with mixed ODEs and algebraic equations.
Computers and Chemical Engineering,

9

(4), 359 – 365.

Cervantes, A., & Biegler, L. T. (1998). Large-scale DAE optimization

using a simultaneous NLP formulation. American Institute of
Chemical Engineering Journal,

44

(5), 1038 – 1050.

Denham, W., & Bryson, A. (1964). Optimal programming problems

with inequality constraints. AIAA Journal,

2

, 25 – 34.

Dixon, L., & Bartholomew Biggs, M. (1981). Adjoint-control trans-

formations for solving practical optimal control problems. Optical
Control Application and Methods,

2

, 365 – 381.

Dreyfus, S. (1962). Variational problems with state variable inequal-

ity constraints. Journal of Mathematical and Analytical Applica-
tions,

4

, 297 – 308.

Duff, I. S. (1981). On algorithm for obtaining a maximum transver-

sal. ACM Transactions on Mathematical Software,

7

(3), 315 – 330.

Feehery, W. F. (1998). Dynamic optimization with path constraints.

Ph.D. thesis, Massachusetts Institute of Technology.

Feehery, W. F., Tolsma, J. E., & Barton, P. I. (1997). Efficient

sensitivity analysis of large-scale differential – algebraic systems.
Applied Numerical Mathematics,

25

(1), 41 – 54.

Gunn, D. J., & Thomas, W. J. (1965). Mass transport and chemical

reaction in multifunctional catalysts. Chemical and Engineering
Science,

20

, 89 – 100.

Ray, W. H. (1981). Ad

6anced process control. New York: McGraw-

Hill.

Haug, E. J., & Ehle, P. E. (1982). Second-order design sensitivity

analysis of mechanical system dynamics. International Journal for
Numerical Methods in Engineering,

18

(1), 1699 – 1717.

Huang, Y. J., Reklaitis, G. V., & Venkatasubramanian, V. (2000).

Dynamic optimization based fault accommodation. Computers
and Chemical Engineering,

24

(2 – 7), 439 – 444 Proceedings of the

process systems engineering conference, PSE 2000, Keystone, CO.

Logsdon, J. S. (1990). Efficient determination of optimal control

profiles for differential algebraic systems. Ph.D. thesis, Carnegie
Mellon University, Pittsburgh, PA.

Maly, T., & Petzold, L. R. (1996). Numerical method and software

for sensitivity analysis of differential – algebraic systems. Applied
Numerical Mathematics,

20

, 57 – 79.

Mattsson, S. E., & So¨derlind, G. (1993). Index reduction in differen-

tial – algebraic equations using dummy derivatives. SIAM Journal
on Scientific Computing,

14

(3), 677 – 692.

Miele, A., Mohanty, B. P., Venkataraman, P., & Kuo, Y. M. (1982).

Numerical solution of minimax problems of optimal control ii.
Journal of Optimization Theory and Applications,

38

(1), 111 – 135.

Morison, K. R. (1984). Optimal control of processes described by

systems of differential – algebraic equations. Ph.D. thesis, Univer-
sity of London, London, UK.

Neuman, C. P., & Sen, A. (1973). A suboptimal control algorithm for

constrained problems using cubic splines. Automatica,

9

(5), 601 –

613.

Rabitz, H., Kramer, M., & Dacol, D. (1983). Sensitivity analysis in

chemical-kinetics. Annual Re

6iew Physical Chemistry,

34

, 419 –

461.

Y.J. Huang et al.

/

Computers and Chemical Engineering

26 (2002) 863 – 873

873

San, K.-Y., & Stephanopoulos, G. (1989). Optimization of fed-batch

penicillin fermentation: a case of singular optimal control with
state constraints. Biotechnology and Bioengineering,

34

, 72 – 78.

Tanartkit, P., & Biegler, L. T. (1995). Stable decomposition for

dynamic optimization. Industrial Engineering and Chemistry Re-
search,

34

(4), 1253 – 1266.

Tanartkit, P., & Biegler, L. T. (1997). A nested, simultaneous ap-

proach for dynamic optimization problems — II: the outer prob-

lem. Computers and Chemical Engineering,

21

(12), 1365 – 1388.

Vasantharajan, S., & Biegler, L. T. (1990). Simultaneous strategies for

optimization of differential – algebraic systems with enforcement
of error criteria. Computers and Chemical Engineering,

14

(10),

1083 – 1100.

Vassiliadis, V. (1993). Computational solution of dynamic optimization

problems with general differential – algebraic constraints. Ph.D. the-
sis, University of London, London, UK.