- English
- français

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

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

**ẋ**(

*t*) =

**f**(

**x**(

*t*),

**u**(

*t*)),

**x**(

*t*

_{0}) =

**x**

_{0}, the cost function

*φ*(

*t*,

_{f}**x**(

*t*)), the

_{f}*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*, or one of its time derivatives, to the states, the inputs or the time derivatives of the inputs, thus resulting in an

_{j}*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,…,

*ρ*-1,

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

**x**

^{u}*being the*

^{ⱼ}*ρ*-dimensional vector of states that can be influenced by manipulating

_{j}*u*, and

_{j}**f**

^{u}*(*

^{ⱼ}**x**,

**u**) such that

**ẋ**

^{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*when it is not determined by an active path constraint, the function det(ℳ

_{j}

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

_{ⱼ}*u*or one of its time derivatives appears in d

_{j}

^{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*) +

**Ť**

_{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,

**Ť**

*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

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

** 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

*ẋ*=

_{j}*f*(

_{j}**x**,

**u**) :=

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

_{j}*ρ*:=

_{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*and its time derivatives, thus resulting in a finite number of solutions, and typically a single solution that satisfies the condition in Case 3.

_{j}Let us define the state vector **x̌*** _{j}* as the complement of the state

*x*(all states

_{j}**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 determined by an active path constraint:

_{j} 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*or one of its time derivatives when it is not determined by an active path constraint.

_{j}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, *Comp. Chem. Eng.*, 27:1-26, 2003.

[2] D. Rodrigues, S. Srinivasan, J. Billeter, and D. Bonvin. Variant and invariant states for chemical reaction systems, *Comp. Chem. Eng.*, 73:23-33, 2015.

[3] B. Chachuat, B. Srinivasan, and D. Bonvin. Adaptation strategies for real-time optimization, *Comp. Chem. Eng.*, 33:1557-1567, 2009.

- EPFL-CONF-228247

#### Reference

Record created on 2017-05-16, modified on 2017-05-17