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

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

dominique.bonvin@epfl.ch

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_x states x(t) described by the differential
equations ẋ(t) = f(x(t),u(t)), x(t₀)
= x₀, the cost function φ(t_f,x(t_f)),
the n_t terminal constraints ψ(t_f,x(t_f))
≤ 0, the n_g mixed path constraints g(x(t),u(t))
≤ 0 and the n_h first-order pure-state constraints 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_j is obtained such that ∂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⁽ⁿ⁾u⁽ⁿ⁺¹⁾,

  Δ_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 ρ_j-dimensional
vector of states that can be influenced by manipulating u_j,
and 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_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ⱼ))/dt^rⱼ, for some r_j. Let u_j^(ξⱼ) be the
highest-order time derivative of u_j that appears in d^rⱼ(det(ℳ_uⱼ))/dt^rⱼ. Then, the optimal input u_j^(ξⱼ) is obtained such
that d^rⱼ(det(ℳ_uⱼ))/dt^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^Tx_r(t)
+ W_in x_in(t) + n₀,

  Q(t) = (-ΔH)^Tx_r(t)
+ Ť_in^T x_in(t)
+ x_ex(t) + Q₀,

where n₀ is the S-dimensional
vector of initial numbers of moles, Q₀ 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, Ť_in is the 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_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 ẋ_j
= f_j(x,u) := u_j, it is possible to define the following vectors of
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_j), and the vector f̌_j(x,u)
as the corresponding complement of f_j(x,u).
Then, one can prove that, when the optimal input u_j is not
determined by an active path constraint:

  1. For reactors with a single independent reaction, the
input u_j is determined by

  d(det(ℳ_uâ±¼))/dt = ∂(∂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.