Infoscience

Conference paper

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 nu piecewise-continuous inputs u(t), nx states x(t) described by the differential equations (t) = f(x(t),u(t)), x(t0) = x0, the cost function φ(tf,x(tf)), the nt terminal constraints ψ(tf,x(tf)) ≤ 0, the ng mixed path constraints g(x(t),u(t)) ≤ 0 and the nh 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 uj be one element of u. The goal is to find an expression that relates the optimal input uj, 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 uj is determined by the active path constraint gk(x,u) = 0.

  2. The optimal input uj is determined by the active path constraint hk(x) ≤ 0, and uj is obtained such that ∂hk/∂x(x) f(x,u) = 0.

  3. Otherwise, the optimal input uj is determined by the condition det(ℳu) = 0, where

  ℳu := [∂fu/∂uj(x,u)  Δjfu/∂uj(x,u)  ⋯  Δj ρ-1fu/∂uj(x,u)],

and the operators Δj,…, Δj ρ-1 are defined as

  Δj v := ∂v/∂x f(x,u) - ∂fu/∂xu(x,u) v + n≥0v/∂u(n) u(n+1),

  Δj l v := Δjj l-1 v),  if l = 2,…, ρj-1,

for any vector field v of dimension ρj, with xu being the ρj-dimensional vector of states that can be influenced by manipulating uj, and fu(x,u) such that u = fu(x,u).

However, the input uj and its time derivatives may not appear explicitly in the function det(ℳu). Hence, as a general approach to find the optimal input uj when it is not determined by an active path constraint, the function det(ℳu) is subject to time differentiation until uj or one of its time derivatives appears in dr(det(ℳu))/dtr, for some rj. Let uj(ξ) be the highest-order time derivative of uj that appears in dr(det(ℳu))/dtr. Then, the optimal input uj(ξ) is obtained such that dr(det(ℳu))/dtr = 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 uin(t) is the p-dimensional vector of inlet flowrates, and qex(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) = NT xr(t) + Win xin(t) + n0,

  Q(t) = (-ΔH)T xr(t) + ŤinT xin(t) + xex(t) + Q0,

where n0 is the S-dimensional vector of initial numbers of moles, Q0 is the initial heat, N is the R×S stoichiometric matrix, ΔH is the R-dimensional vector of heats of reaction, Win is the S×p inlet-composition matrix, Ťin is the p-dimensional vector of inlet specific enthalpies, xr(t) is the R-dimensional vector of extents of reaction, xin(t) is the p-dimensional vector of extents of inlet, and xex(t) is the extent of heat exchange.

The state vector of dimension nx := R + p + 1 is

  x(t) := [xr(t)Txin(t)Txex(t)]T,

while the input vector of dimension nu := p + 1 is

  u(t) := [uin(t)Tqex(t)]T.

The dynamic equations can be written compactly as

  (t) = f(x(t),u(t)), x(t0) = 0R+p+1,

  f(x(t),u(t)) := [rv(x(t))Tu(t)T]T,

where rv(x(t)) is the R-dimensional vector of reaction rates. In batch and semi-batch reactors, with j = fj(x,u) := uj, it is possible to define the following vectors of dimension ρj := R+1:

  xu := [xrTxj]T,

  fu(x,u) := [rv(x)Tuj]T.

Note that, since this system is input-affine, det(ℳu) and its time derivatives are polynomial functions of uj 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 j as the complement of the state xj (all states x except xj), and the vector j(x,u) as the corresponding complement of fj(x,u). Then, one can prove that, when the optimal input uj is not determined by an active path constraint:

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

  d(det(ℳu))/dt = ∂(∂rv/∂xj(x))/∂xj uj + ∂(∂rv/∂xj(x))/∂j j(x,u),

since uj and its time derivatives do not appear in

  det(ℳu) = ∂rv/∂xj(x).

  2. For reactors with 2 independent reactions, the input uj is determined by

  det(ℳu) = det([∂rv/∂xj(x)  ∂(∂rv/∂xj(x))/∂xj]) uj

   + det([∂rv/∂xj(x)  -∂rv/∂xr(x) ∂rv/∂xj(x) + ∂(∂rv/∂xj(x))/∂j j(x,u)])

One can use symbolic computation software to evaluate the function det(ℳu) and its time derivatives and obtain the optimal input uj 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, 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.

    Reference

    • EPFL-CONF-228247

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

Related material