# (430b) Optimal Control Laws for Batch and Semi-Batch Reactors Using the Concept of Extents

#### AIChE Annual Meeting

#### 2017

#### 2017 Annual Meeting

#### Computing and Systems Technology Division

#### In Honor of Christos Georgakis' 70th Birthday

#### Tuesday, October 31, 2017 - 3:33pm to 3:51pm

**Optimal Control Laws for Batch and Semi-batch Reactors
Using the Concept of Extents**

**D. Rodrigues, J. Billeter and D. Bonvin**

Laboratoire dâ€™Automatique, Ecole Polytechnique FÃ©dÃ©rale de Lausanne

EPFL â€“ Station 9, 1015 Lausanne, Switzerland

In the context of dynamic optimization of batch and

semi-batch reactors, this contribution presents a method that uses the concept

of extents to generate solutions that satisfy the necessary conditions of optimality

given by Pontryaginâ€™s maximum principle. The method is

divided in two parts. In the first part, the reactor model is written in terms

of decoupled extents, and adjoint-free optimal control laws are generated for

all possible types of arcs that may occur in the optimal solution. In the

second part, the correct sequence of arcs is determined, and, for each sequence,

the optimal switching times and initial conditions are computed numerically.

The optimal control problems (OCPs) are formulated in Mayer

form, with *n _{u}* piecewise-continuous inputs

**u**(

*t*),

*n*states

_{x}**x**(

*t*) described by the differential

equations

**xÌ‡**(

*t*) =

**f**(

**x**(

*t*),

**u**(

*t*)),

**x**(

*t*

_{0})

=

**x**

_{0}, the cost function

*Ï†*(

*t*,

_{f}**x**(

*t*)),

_{f}the

*n*terminal constraints

_{t}**Ïˆ**(

*t*,

_{f}**x**(

*t*))

_{f}â‰¤

**0**, the

*n*mixed path constraints

_{g}**g**(

**x**(

*t*),

**u**(

*t*))

â‰¤

**0**and the

*n*first-order pure-state constraints

_{h}**h**(

**x**(

*t*))

â‰¤

**0**. Note that this class of OCPs is not restrictive in most dynamic optimization

problems dealing with reactors.

The optimal input trajectories are composed of a (typically finite) number of arcs.

For each arc and for each input, the optimal input is determined by either an

active path constraint or a condition that expresses a physical compromise that

depends exclusively on the dynamics of the system [1].

Let the input *u _{j}* be one element of

**u**.

The goal is to find an expression that relates the optimal input

*u*,

_{j}or one of its time derivatives, to the states, the inputs or the time

derivatives of the inputs, thus resulting in an

*adjoint-free*optimal control law. For each arc, one of the following

cases is possible:

1. The optimal

input *u _{j}* is determined by the active path constraint

*g*(

_{k}**x**,

**u**)

= 0.

2. The

optimal input *u _{j}* is determined by the active path constraint

*h*(

_{k}**x**)

â‰¤ 0, and

*u*is obtained such that âˆ‚

_{j}*h*/âˆ‚

_{k}**x**(

**x**)

**f**(

**x**,

**u**)

= 0.

3. Otherwise,

the optimal input *u _{j}* is determined by the condition det(â„³

_{u}*) = 0, where*

_{â±¼} â„³_{u}* _{â±¼}*Â := [âˆ‚

**f**

^{u}*/âˆ‚*

^{â±¼}*u*(

_{j}**x**,

**u**)

Î”

*âˆ‚*

_{j}**f**

^{u}*/âˆ‚*

^{â±¼}*u*(

_{j}**x**,

**u**)

â‹¯

Î”

_{j}

_{}

^{Ï}

^{â±¼}^{-1}

âˆ‚

**f**

^{u}*/âˆ‚*

^{â±¼}*u*(

_{j}**x**,

**u**)],

and the operators Î”* _{j}*,â€¦, Î”

_{j}

_{}

^{Ï}

^{â±¼}^{-1}

are defined as

Î”_{j}**v** := âˆ‚**v**/âˆ‚**x**

**f**(**x**,**u**) - âˆ‚**f**^{u}* ^{â±¼}*/âˆ‚

**x**

^{u}*(*

^{â±¼}**x**,

**u**)

**v**+

**âˆ‘**

_{n}_{â‰¥0}âˆ‚

**v**/âˆ‚

**u**

^{(n)}

**u**

^{(n+1)},

Î”_{j}_{}^{l}**v** := Î”_{j}

(Î”_{j}_{}^{l}^{-1} **v**),

if *l* = 2,â€¦, *Ï _{j}*-1,

for any vector field **v** of dimension *Ï _{j}*,

with

**x**

^{u}*Â being the*

^{â±¼}*Ï*-dimensional

_{j}vector of states that can be influenced by manipulating

*u*,

_{j}and

**f**

^{u}*(*

^{â±¼}**x**,

**u**) such that

**xÌ‡**

^{u}*Â =*

^{â±¼}**f**

^{u}*(*

^{â±¼}**x**,

**u**).

However, the input *u _{j}* and its time

derivatives may not appear explicitly in the function det(â„³

_{u}*). Hence, as a general approach to find the optimal input*

_{â±¼}*u*

_{j}when it is not determined by an active path constraint, the function det(â„³

_{u}*) is subject to time differentiation until*

_{â±¼}*u*

_{j}or one of its time derivatives appears in d

^{r}*(det(â„³*

^{â±¼}

_{u}*))/d*

_{â±¼}*t*

^{r}*, for some*

^{â±¼}*r*. Let

_{j}*u*

_{j}^{(Î¾}

^{â±¼}^{)}be the

highest-order time derivative of

*u*that appears in d

_{j}

^{r}*(det(â„³*

^{â±¼}

_{u}*))/d*

_{â±¼}*t*

^{r}*. Then, the optimal input*

^{â±¼}*u*

_{j}^{(Î¾}

^{â±¼}^{)}is obtained such

that d

^{r}*(det(â„³*

^{â±¼}

_{u}*))/d*

_{â±¼}*t*

^{r}*Â = 0.*

^{â±¼}The mass and heat balances for batch and semi-batch reactors

can be written using the concept of extents [2]. Let us consider a homogeneous

batch or semi-batch reactor with *R* independent reactions and *p*

independent inlets (*p *= 0 for batch reactors), where **u*** _{in}*(

*t*)

is the

*p*-dimensional vector of inlet flowrates, and

*q*

*(*

_{ex}*t*)

is the exchanged heat power. The numbers of moles

**n**(

*t*) and the

heat

*Q*(

*t*) can be expressed as a linear combination of

extents, according to

** n**(*t*) = **N**^{T}**x*** _{r}*(

*t*)

+

**W**

_{in}**x**

*(*

_{in}*t*) +

**n**

_{0},

* Q*(*t*) = (-Î”**H**)^{T}**x*** _{r}*(

*t*)

+

**TÌŒ**

_{in}^{T}

**x**

*(*

_{in}*t*)

+

*x*(

_{ex}*t*) +

*Q*

_{0},

where **n**_{0} is the *S*-dimensional

vector of initial numbers of moles, *Q*_{0} is the initial heat, **N** is the *R*Ã—*S* stoichiometric matrix, Î”**H** is the *R*-dimensional vector of heats of

reaction, **W*** _{in}* is the

*S*Ã—

*p*inlet-composition matrix,

**TÌŒ**

*is the*

_{in}*p*-dimensional vector of inlet specific

enthalpies,

**x**

*(*

_{r}*t*) is the

*R*-dimensional

vector of extents of reaction,

**x**

*(*

_{in}*t*) is the

*p*-dimensional

vector of extents of inlet, and

*x*(

_{ex}*t*) is the extent

of heat exchange.

The state vector of dimension *n _{x}* :=

*R*+

*p*+ 1 is

** x**(*t*) := [**x*** _{r}*(

*t*)

^{T}

**x**

*(*

_{in}*t*)

^{T}

*x*(

_{ex}*t*)]

^{T},

while the input vector of dimension *n _{u}* :=

*p*+ 1 is

** u**(*t*) := [**u*** _{in}*(

*t*)

^{T}

*q*

*(*

_{ex}*t*)]

^{T}.

The dynamic equations can be

written compactly as

**xÌ‡**(*t*) = **f**(**x**(*t*),**u**(*t*)), **x**(*t*_{0}) = **0**_{R}_{+p+1},

by defining

** f**(**x**(*t*),**u**(*t*)) := [**r*** _{v}*(

**x**(

*t*))

^{T}

**u**(

*t*)

^{T}]

^{T},

where **r*** _{v}*(

**x**(

*t*)) is

the

*R*-dimensional vector of reaction rates. In batch and semi-batch reactors, with

*xÌ‡*

_{j}=

*f*(

_{j}**x**,

**u**) :=

*u*, it is possible to define the following vectors of

_{j}dimension

*Ï*:=

_{j}*R*+1:

** x**^{u}* ^{â±¼}*Â := [

**x**

_{r}^{T}

*x*]

_{j}^{T},

** f**^{u}* ^{â±¼}*(

**x**,

**u**) := [

**r**

*(*

_{v}**x**)

^{T}

*u*]

_{j}^{T}.

Note that, since this system is input-affine, det(â„³_{u}* _{â±¼}*) and its time derivatives are polynomial functions of

*u*

_{j}and its time derivatives, thus resulting in a finite number of solutions, and typically a single solution

that satisfies the condition in Case 3.

Let us define the state vector **xÌŒ**_{j}

as the complement of the state *x _{j}* (all states

**x**except

*x*), and the vector

_{j}**fÌŒ**

*(*

_{j}**x**,

**u**)

as the corresponding complement of

*f*(

_{j}**x**,

**u**).

Then, one can prove that, when the optimal input

*u*is not

_{j}determined by an active path constraint:

1. For reactors with a single independent reaction, the

input *u _{j}* is determined by

d(det(â„³_{u}* _{â±¼}*))/d

*t*= âˆ‚(âˆ‚

**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**))/âˆ‚

*x*

_{j}*u*+ âˆ‚(âˆ‚

_{j}**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**))/âˆ‚

**xÌŒ**

_{j}**fÌŒ**

*(*

_{j}**x**,

**u**),

since *u _{j}* and its time derivatives do not

appear in

det(â„³_{u}* _{â±¼}*) = âˆ‚

**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**).

2. For reactors with 2 independent reactions, the input *u _{j}*

is determined by

det(â„³_{u}* _{â±¼}*) = det([âˆ‚

**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**)

âˆ‚(âˆ‚

**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**))/âˆ‚

*x*])

_{j}*u*

_{j} + det([âˆ‚**r*** _{v}*/âˆ‚

*x*(

_{j}**x**)

-âˆ‚

**r**

*/âˆ‚*

_{v}**x**

*(*

_{r}**x**)

âˆ‚

**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**) + âˆ‚(âˆ‚

**r**

*/âˆ‚*

_{v}*x*(

_{j}**x**))/âˆ‚

**xÌŒ**

_{j}**fÌŒ**

*(*

_{j}**x**,

**u**)])

One can use

symbolic computation software to evaluate the function det(â„³_{u}* _{â±¼}*) and its time derivatives and obtain the optimal input

*u*

_{j}or one of its time derivatives when it is not determined by an active path

constraint.

The advantage of the proposed approach is that it reduces

the set of possible arcs to a finite

number of possibilities. This, in turn, results in a finite number of arc

sequences if one assumes an upper bound on the number of arcs present in the

optimal solution. Hence, instead of solving the original infinite-dimensional

problem, one can simply perform numerical optimization for each arc sequence,

using the switching times between arcs and the initial conditions as decision variables.

Note that the dynamic state equations

can be integrated forward in time, since it is possible to evaluate the

corresponding inputs without knowledge of the adjoint variables. Once the

forward integration is complete, one can integrate backward in time to obtain

the corresponding adjoint variables, which enables the computation of the

gradients with respect to the switching times and initial conditions of the

arcs.

The proposed approach will be

illustrated to maximize the final quantity

of product in an acetoacetylation reaction, subject to constraints on the final

concentration of reactants and by-products [3].

[1] B. Srinivasan, S. Palanki, and D. Bonvin. Dynamic

optimization of batch processes: I. Characterization of the nominal solution, *Comput.
Chem. Eng.*, 27:1-26, 2003.

[2] D. Rodrigues, S. Srinivasan, J. Billeter, and D.

Bonvin. Variant and invariant states for chemical reaction systems, *Comput.
Chem. Eng.*, 73:23-33, 2015.

[3] B. Chachuat, B. Srinivasan, and D.

Bonvin. Adaptation strategies for real-time optimization, *Comput. Chem. Eng.*,

33:1557-1567, 2009.