Robust optimization of stiff delayed systems: application to a fluid catalytic cracking unit

Abstract

We apply a robust steady state optimization method for stiff delay differential equations to the economic optimization of a fluidized catalytic cracking unit. Stiff systems of differential equations appear in this case due to the different time scales in the gas and fluid phase. Delays result from the catalyst hold-ups in the standpipes connecting riser and regenerator. We show that the proposed robust optimization method can cope with stiffness and delays. Moreover, the proposed method is capable of simultaneously optimizing the process parameters and tuning controller parameters.

1 Introduction

Fluidized catalytic cracking (FCC) increases the gasoline yield and therefore it is one of the key technologies in crude oil processing. The process of fluidized catalytic cracking breaks up the fractions of heavy hydrocarbons into lighter and more valuable molecules. The cracking process is dominated by a catalytic reaction taking place at the surface of a solid catalyst which is used in a powder form to achieve a large active surface. If the cracking reaction takes place under plug flow conditions, the contact time of feedstock and catalyst can be controlled well. Excess carbon accumulates on the catalyst surface as coke, deactivating it by time. Used catalyst is led into the regenerator, where air provides oxygen to burn off the coke and keeps the catalyst fluidized. The regenerated catalyst is then ready for use in the cracking process.

FCC process dynamics can exhibit state multiplicity (Arbel et al. 1995; Pinheiro et al. 2012), which is attended by the existence of fold bifurcations (Elnashaie et al. 1995). Hopf bifurcations also occur in FCC under certain conditions (Alhumaizi and Elnashaie 1997; Elnashaie et al. 1995). The presence of bifurcations complicates steady-state optimizations because an apparently optimal steady state may turn out to lie on an unstable solution branch of the steady-state manifold. Such an optimization result cannot be implemented on a real plant without further measures for stabilization. There exist, however, optimization methods that enforce stability at the optimal steady state.

Stability of a steady state can be achieved in an optimization by penalizing the largest eigenvalue real part (Vanbiervliet et al. 2008) or its expected value (Fenzi and Michiels 2017). In combination with an economic profit function to be maximized, such a treatment of stability will result in a multi-objective optimization, calling for Pareto optimization or a weighting of the optimization objectives.

Alternatively, stability requirements can be enforced with constraints instead of with the cost function. So-called normal vector methods have been successfully applied to various optimization problem types and processes, e. g., to the reaction section of a HDA process (Mönnigmann and Marquardt 2005), a CSTR with recycling delay (Kastsian and Mönnigmann 2013) and to finding optimal and robustly stable limit cycles of chemical reactions (Kastsian and Mönnigmann 2014). Applications aiming for the robust satisfaction of other criteria like robustly limiting overshooting (Gerhard et al. 2008) or robust synchronization of subsystems are also possible (Otten et al. 2018). Normal vector methods can also be used for an automated controller tuning of closed-loop nonlinear systems (Hahn et al. 2008).

The following two aspects render the application of normal vector methods to fluid catalyzed cracking challenging: Firstly, in a fluidized cracking unit, the standpipes connecting regenerator and riser have to hold up a sufficient amount of catalyst to seal off the oxygenated atmosphere of the regenerator and the oxygen-free atmosphere within the riser. This hold-up is known to cause a transport delay which is often neglected in dynamical models (Han and Chung 2001a) for simplicity. Secondly, the typical FCC process model is based on the assumption of fast gas dynamics which justifies the use of algebraic equations instead of differential equations. This results in differential algebraic equations (DAEs) (Ali et al. 1997). If the gas dynamics are instead modeled with differential equations, the different time scales result in a system of stiff differential equations. If delays must be taken into account, a system of delay differential algebraic equations or a system of stiff delay differential equations must be treated. The numerical treatment of these classes is still subject to research. Specialized numerical integrators, which are necessary to run simulations of the FCC model treated here, have been developed only recently (DelayDiffEq.jl; Rackauckas and Nie 2017).

This paper is structured as follows. Section 2 introduces a dynamical model of a fluidized catalytic cracker. A system of delay differential equations results that are stiff due to different time scales of gas phase and solid phase processes. The steady-state optimization problem is presented in Sect. 3 to motivate constraints guaranteeing robust stability. A suitable method to achieve robust stability is introduced in Sect. 4 and applied to find an optimal steady state of the FCC. A summary and conclusions are stated in Sect. 6.

2 Modeling

The model is based on a model proposed by Ali and Rohani (1997). It focuses on the components in which reactions take place, i. e., the riser and the regenerator. Their properties are sketched in the following sections.

2.1 Riser

Plug fow is beneficial for the cracking reaction (Shaikh et al. 2008). Modeling a plug flow results in a partial differential equation. The classical approach to analyzing FCC dynamics is to assume the plug flow reactor dynamics fast enough compared to the regenerator dynamics to justify a steady state assumption (Ali and Rohani 1997; Ali et al. 1997; Arbel et al. 1995; Shaikh et al. 2008). This eliminates the time derivatives and reduces the partial differential equation for the plug flow reactor to a spatial ordinary differential equation. Alternatively, the underlying PDE can be discretized spatially (Secchi et al. 2001), thus retaining the plug flow reactor dynamics. The latter approach is used here, where a uniform mesh of $L=15$ discrete points is used, which corresponds to approximating the plug flow reactor by a sequence of 15 CSTRs.

The reactions taking place in the riser (and therefore in each CSTR) are modeled by four-lump models, in which the different types of hydrocarbons are divided into four classes: the gas oil A is converted into gasoline B, coke C and light hydrocarbons D. Gasoline is the desired product, coke and light hydrocarbons are byproducts. Gasoline is converted to light hydrocarbons and coke in an undesired side reaction. The differential equations at the discrete point i read

$$\begin{aligned} \dot{y}_{A,l}&= -\frac{F_{gR}}{A_{\text {Ris}}\,\epsilon _{gR}\,\rho _{gR}\,U_{\text {Ris}}}\frac{1}{h}(y_{A,l}-y_{A,l-1}) - \frac{\Phi _l L_{\text {Ris}}}{U_{\text {Ris}}}(k_{AB,l} + k_{AC,l} + k_{AD,l})\,y_{A,l}^2\,, \end{aligned}$$

(1a)

$$\begin{aligned} \dot{y}_{B,l}&= -\frac{F_{gR}}{A_{\text {Ris}}\,\epsilon _{gR}\,\rho _{gR}\,U_{\text {Ris}}}\frac{1}{h}(y_{B,l}-y_{B,l-1}) - \frac{\Phi _l L_{\text {Ris}}}{U_{\text {Ris}}}(k_{BC,l}\,y_{B,l} - k_{AB,l}\, y_{A,l}^2)\,,\end{aligned}$$

(1b)

$$\begin{aligned} \dot{y}_{C,l}&= -\frac{F_{gR}}{A_{\text {Ris}}\,\epsilon _{gR}\,\rho _{gR}\,U_{\text {Ris}}}\frac{1}{h}(y_{C,l}-y_{C,l-1}) + \frac{\Phi _l L_{\text {Ris}}}{U_{\text {Ris}}}(k_{BC,l}\,y_{B,l} + k_{AC,l}\, y_{A,l}^2)\,,\end{aligned}$$

(1c)

$$\begin{aligned} \dot{y}_{D,l}&= -\frac{F_{gR}}{A_{\text {Ris}}\,\epsilon _{gR}\,\rho _{gR}\,U_{\text {Ris}}}\frac{1}{h}(y_{D,l}-y_{D,l-1}) + \frac{\Phi _l L_{\text {Ris}}}{U_{\text {Ris}}}(k_{BD,l}\,y_{B,l} + k_{AD,l}\, y_{A,l}^2)\,,\end{aligned}$$

(1d)

$$\begin{aligned} \dot{T}_{\text {Ris},l}&= -\frac{F_{s}\,c_{ps}+F_{gR}\,c_{prG}}{A_{\text {Ris}}(1-\epsilon _{gR})\,\rho _{gR}\,U_{\text {Ris}}}\frac{1}{h}(T_{\text {Ris},l}-T_{\text {Ris},l-1}) - \frac{c}{h}(T_{\text {Ris},l} - T_\text {out})\nonumber \\&\quad +\frac{\Phi _l \epsilon _{gR}\,\rho _{gR}\,L_{\text {Ris}}}{(1-\epsilon _{gR})\rho _s\, c_{ps} \,U_{\text {Ris}}t_\text {ref}}((\Delta H_{AB}k_{AB,l}+\Delta H_{AC} k_{AC,l}+\Delta H_{AD} k_{AD,l})\,y_{A,l}^2 \nonumber \\&\quad + (\Delta H_{BC}k_{BC,l}+\Delta H_{BD}k_{BD,l})\,y_{B,l})\,,\end{aligned}$$

(1e)

$$\begin{aligned} {{\dot{\Phi }}}_l&= -\frac{F_{gR}}{A_{\text {Ris}}\,\epsilon _{gR}\,\rho _{gR}\,U_{\text {Ris}}}\frac{1}{h}(\Phi _l - \Phi _{l-1}) -\alpha _l \Phi _l\,, \end{aligned}$$

(1f)

where $y_{A_l}$, $y_{B_l}$, $y_{C_l}$ and $y_{D_l}$ describe the concentrations of each lump at the location l scaled to the gas oil concentration at the riser bottom ($y_{A_0}$). $T_{\text {Ris},l}$ is the temperature at this point and $\Phi _l$ is the catalyst activity scaled to the activity at the riser bottom ($\Phi _{0}$). All constants and $y_{A,0}$, $y_{B,0}$, $y_{C,0}$, $y_{D,0}$, $T_{\text {Ris},0}$ and $\Phi _{0}$, which are necessary to consider the boundary condition at $l=1$, are listed in Tables 5, 6 and 7 in “Appendix A”.

At each discrete point, $n_{x,\text {Ris},l} = 6$ states are taken into account, which results in $n_{x,\text {Ris}} = 90$ states to describe the whole riser.

2.2 Regenerator

It is assumed in the regenerator model that nearly all of the catalyst is located at the lower part of the vessel (Ali and Rohani 1997). Three phases must be considered: A solid phase of catalyst particles, a gaseous emulsion of the gases surrounding the particles which keep the catalyst fluidized, and a gaseous bubble phase, consisting of excess gas traveling upwards through the catalyst bed.

The solid phase and the bubble phase can only interact via the emulsion phase. The emulsion phase is well-mixed due to the turbulent character of fluidization. Consequently, a CSTR model suffices.

The interaction of the bubble and emulsion phase is very fast compared to the interaction of solid phase and emulsion phase (Ali and Rohani 1997). The bubble phase is therefore assumed to be at equilibrium with the emulsion phase at all times. Hence, the plug-flow dynamics in the bubble phase can be explicitly solved and substituted into the emulsion phase equations. The resulting differential equations read

$$\begin{aligned} \dot{W}_{cg}&= F_{s}( W_{cr} - W_{cg} ) - \frac{\rho _b A_G L_G }{M_{SG}}( 1 - \epsilon _{bG} ) k_c W_{cg} y_{O_2}C_{O_2f} \,, \end{aligned}$$

(2a)

$$\begin{aligned} \dot{y}_{O_2}&= A_G \left( U_\text {mf} + ( U_{ Reg} - U_\text {mf} )\frac{ 1 - \exp ( -\alpha _I )}{( L_G A_G \epsilon _{dG} )}\right) ( y_{O_2f} - y_{O_2} )\nonumber \\&\quad - \frac{( 1 - \epsilon _{dG} ) \rho _b}{ \epsilon _{dG} }\left( 0.5 k_{COd}^\prime \,y_{CO}\, y_{O_2}^{0.5} + \frac{k_{C}^{\prime \prime \prime }}{M_{Wc}}W_{cg} y_{O_2} \right) \,,\end{aligned}$$

(2b)

$$\begin{aligned} \dot{y}_{CO}&= A_G\left( U_\text {mf} + ( U_{\text {Reg}} - U_\text {mf} )\frac{ 1 - \exp ( -\alpha _I )}{ L_G A_G \epsilon _{dG} }\right) ( y_{COf} - y_{CO} ) \nonumber \\&\quad - \frac{( 1 - \epsilon _{dG} )\rho _b}{ \epsilon _{dG} }\left( k_{COd}^\prime y_{CO}\, y_{O_2}^{0.5} - \frac{ k_C^\prime W_{cg}y_{O_2}}{M_{Wc}} \right) \,,\end{aligned}$$

(2c)

$$\begin{aligned} \dot{y}_{CO_2}&= A_G\left( U_\text {mf} + ( U_{\text {Reg}} - U_\text {mf} )\frac{ 1 - \exp ( -\alpha _I )}{ L_G A_G\epsilon _{dG} }\right) ( y_{CO_2f} - y_{CO_2} ) \nonumber \\&\quad + \frac{ (1 - \epsilon _{dG} )\rho _b }{\epsilon _{dG}} \left( 0.5 k_{COd}^\prime y_{CO} y_{O_2}^{0.5} + \frac{k_{C}^{\prime \prime }}{M_{Wc}}W_{cg} y_{O_2} )\right) \,,\end{aligned}$$

(2d)

$$\begin{aligned} \dot{T}_\text {Reg}&= A_G\left( U_\text {mf} + ( U_{\text {Reg}} - U_\text {mf} )\frac{ (\exp ( -\alpha _H ) - 1 )\rho _{gG}c_{pgG} )}{ M_{SG}c_{ps}} \right) ( T_\text {air} - T_\text {Reg} ) \nonumber \\&\quad + \frac{ F_{s} }{ M_{SG} c_{ps} }( T_{\text {Top}} - T_\text {Reg} ) -c ( T_\text {Reg}-T_\text {out} )\nonumber \\&\quad + \frac{ \rho _b A_G L_G ( 1 - \epsilon _{bG} )C_{O_2f} }{ M_{SG} c_{ps} t_\text {ref}}\left( \frac{ k_c W_{cg} y_{O_2} }{ M_{Wc} }\Delta H_{Rc}+k_{cod}y_{CO} y_{O_2}^{0.5}\Delta H_{Rcod} \right) \,, \end{aligned}$$

(2e)

where $W_{cg}$ describes the coke mass per catalyst mass, $T_{Reg}$ describes the regenerator temperature and $y_{O_2}$, $y_{CO}$ and $y_{CO_2}$ describe the concentration of oxygen, carbon monoxide and carbon dioxide scaled to the oxygen concentration in fresh air feed, respectively.

This submodel contributes $n_{x,\text {Reg}} = 5$ states to the FCC model (cf. Table 1).

2.2.1 Delays

Standpipes filled with catalyst separate the oxygen-rich atmosphere in the regenerator from the flammable hot hydrocarbon atmosphere in the riser. This inflicts transport delays on the dynamics, as the heat energy and coke-bearing catalyst have to pass a standpipe when transferred from the riser to the regenerator ($\tau _1$) and vice versa ($\tau _2$). This transport delay lies in the magnitude of seconds and depends on the geometry of the FCC (Han and Chung 2001a).

The numerical values of the delays given in Table 7 correspond to assumed standpipe lengths of $8 \text { to } 10\,{\mathrm {m}}$ and a diameter of $0.14\,{\mathrm {m}}$. The coke-bearing catalyst from the riser top enters the regenerator with a delay $\tau _1$ which affects $W_{cr}$ and $T_\text {Top}$. The regenerated catalyst enters the riser with a delay $\tau _2$, which affects $T_{\text {Ris},0}$.

2.3 Control structure

Controller parameters are tuned by the normal vector based steady-state optimization. The controller structure must be chosen in advance, however. Alhumaizi and Elnashaie (1997) compared six control structures, of which the following is suited best to influence the occurrence of bifurcations and it offers flexibility in the handling of different feedstock (Arbel et al. 1997). The air feed temperature $T_{\text {air}}$ is manipulated to control the regenerator temperature, while the catalyst circulation $F_s$ rate is manipulated to control the riser temperature,

$$\begin{aligned} T_{\text {air}}&= K_\text {Reg}( {\hat{T}}_\text {Reg,SP}- T_\text {Reg}) + T_{\text {air},0}\,, \end{aligned}$$

(3a)

$$\begin{aligned} F_{\text {s}}&= K_\text {Ris}( {\hat{T}}_\text {Ris,SP} - T_\text {Ris}) + F_{\text {s},0}. \end{aligned}$$

(3b)

For nonzero controller gains, the offsets $T_{\text {air},0}$ and $F_{\text {s},0}$ can be omitted by substituting ${\hat{T}}_\text {reg,SP} = T_\text {reg,SP} - \frac{T_{\text {air},0}}{K_\text {Reg}}$ and ${\hat{T}}_\text {Ris,SP} = T_\text {Ris,SP} - \frac{F_{\text {s},0}}{K_\text {Ris}}$, respectively.

2.3.1 Combined model

The process is shown in Fig. 1. It includes the delays and controllers described above.

Fig. 1

States and parameters of the FCC with its controllers. The parameters shown are optimization variables. Parameters not listed here are constant during the optimization. They are given in “Appendix A”

We denote the delay differential equations that result from combining (1), (2) and (3) as

$$\begin{aligned} \dot{x}(t)=f(x(t),x(t-\tau _1),x(t-\tau _2),\alpha ,p)\, \end{aligned}$$

(4)

for brevity. Because of the delays discussed in Sect. 2.2.1 these equations not only depend on the state variables x(t) at the current time t, but also on $x(t- \tau _1)$ and $x(t-\tau _2)$. Assume, for example, the catalyst requires a time $\tau _1$ to pass the standpipe when being transported from the riser to the regenerator. The catalyst feed temperature at the regenerator ($T_{Top}$, see (2e)) then only assumes the riser exit temperature $T_{Ris, L}$ after this delay $\tau _1$. This results in the relation $T_{Top}(t)= T_{Ris, L}(t- \tau _1)$, which implies both the state variable $T_{Ris, L}(t)$ and $T_{Ris, L}(t- \tau _1)$ appear in (4). Similarly, the delay $\tau _2$ explained in Sect. 2.2.1 introduces a dependency on $x(t- \tau _2)$ in (4). The system of delay differential equations (4) consists of $n_x=95$ state variables x(t) and has seven parameters, which are listed in Tables 1 and 2, respectively. Parameters that are and are not subject to uncertainty are denoted by $\alpha$ and p, respectively. We assume all but the controller parameters to be subject to uncertainty ($\alpha \in {\mathbb {R}}^{n_\alpha }$, $n_\alpha = 5$, $p\in {\mathbb {R}}^{n_p}$, $n_p= 2$, cf. Table 2).

The uncertainties are assumed to originate from measurement uncertainties. However, the uncertainties of the temperature setpoints require an explanation. They represent the uncertainties of the temperature measurements,

$$\begin{aligned} K_i ( (T_i \pm \Delta T_i) - T_{i,SP}) = K_i ( T_i - (T_{i,SP}\mp \Delta T_i)). \end{aligned}$$

Table 1 FCC model states

Table 2 FCC optimization parameters

An initial analysis of this model reveals that, for some parameter combinations, multiple steady states exist. This state multiplicity is visualized in Fig. 2. This renders the optimization more difficult, as the state multiplicity also implies the presence of unstable steady states.

Fig. 2

Three steady states exist for typical gas oil feed rates (e. g., $F_{gR} = 19.95\, {\mathrm {kg}}/{\mathrm {s}}$ at $T_{gR} = 494\,{\mathrm {K}}$). When multiplicity occurs, steady states with the lowest and highest temperatures are stable (solid lines), while the middle one is unstable (dashed lines). The remaining parameters are $F_\text {air}=16\, {\mathrm {kg}}/{\mathrm {s}}$, $T_\text {Reg,SP} = 1000\,{\mathrm {K}}$, $T_\text {Ris,SP} = 850\,{\mathrm {K}}$, $K_\text {Reg} = 1$ and $K_\text {Ris} = 1$. a Regenerator temperature T_Reg at steady states for various gas oil feeds F_gR. The isola causing multiple steady states only exists for a range of gas oil feeds. b Regenerator temperature T_Reg at steady states for various gas oil feed temperatures T_gR. Gas oil feed temperature axis extends to beyond usual process conditions; shown for completeness only

3 Optimization problem statement

We optimize the FCC with respect to the profit function

$$\begin{aligned} \phi = (P_B\,y_{B,L} + P_D\,y_{D,L}-P_A (1-y_{A,L}))F_{gR}\,, \end{aligned}$$

(5)

where the prices $P_A$, $P_B$ and $P_D$ of gas oil, gasoline and light hydrocarbons, respectively, are given in Table 3. The quantities $y_{A,L}$, $y_{B,L}$ and $y_{D,L}$ are the gas oil concentration, gasoline concentration and light hydrocarbon concentration, respectively, at the riser top. $F_{gR}$ is, as stated in Table 2, the gas oil feed.

Table 3 Given prices for gas oil and gasoline are based on stock market prices, converted to $\$/{\mathrm {ton}}$ and rounded to $10\$$. The price for light hydrocarbons is chosen to lie between prices of natural gas and gas oil

Combining the cost function (5) with constraints that enforce a steady state results in the optimization problem

$$\begin{aligned}&\max _{x,\alpha ,p}\,\phi \end{aligned}$$

(6a)

$$\begin{aligned}&\text {s.\,t. } \,f(x,x,x,\alpha ,p) = 0 \end{aligned}$$

(6b)

$$\begin{aligned}&\alpha \,\text {in feasible intervals stated in Table~2} \end{aligned}$$

(6c)

$$\begin{aligned}&2\,{\mathrm {kg}}/{\mathrm {s}}\le F_s\le 200\,{\mathrm {kg}}/{\mathrm {s}} \end{aligned}$$

(6d)

$$\begin{aligned}&300\,{\mathrm {K}}\,\le T_\text {air}\,\le 2500\,{\mathrm {K}}\,, \end{aligned}$$

(6e)

where state variables are as given in Table 1, $\alpha = [F_\text {air}, F_{gR}, T_{gR},$ $T_\text {Reg,SP},T_\text {Ris,SP} ]^\prime$ and $p=[K_\text {Reg},K_\text {Ris}]^\prime$. The inequalities (6d) and (6e) are simple bounds on the catalyst circulation rate $F_s$ and the air feed temperature $T_\text {air}$. By substituting (3), $F_s$ and $T_\text {air}$ can be expressed in terms of the remaining variables and parameters.

The solution to this optimization may be unstable, as the stability of the steady state is not taken into account in (6). Without any further stabilization effort, such an unstable steady state is of little practical use. The cost function value that results for (6) will, however, serve as a benchmark for the optimization guaranteeing robust stability (see Table 4 in Sect. 5). Constraints that ensure robust stability are introduced in the next section.

4 Normal vector constraints for robust stability

We introduce the central idea of the normal vector constraints informally first. Essentially, these constraints enforce a safe distance to any critical boundary on the steady state manifold of the optimized system. The term “critical boundary” may refer to any boundary in the space of the parameters $\alpha$ that separates steady states with desired dynamical properties from those with undesired dynamical properties (Mönnigmann and Marquardt 2002). In the present paper, we are interested in boundaries that separate stable from unstable steady states. These boundaries consist of fold bifurcation points here.

Fig. 3

Sketch of a steady state manifold with a critical boundary defined by $G(\cdot )=0$. The projection of the critical manifold locally separates regions with desired from those with undesired behavior (such as stable from unstable steady states)

Figure 3 sketches a manifold of steady states, i.e., any point on the surface represents a pair $x^{(0)}$, $\alpha ^{(0)}$ such that $f(x^{(0)}, \dots , x^{(0)}, \alpha ^{(0)}, p)= 0$ for arbitrary but fixed parameter values p. The stability properties of such a steady state can be characterized with the eigenvalues of the Jacobian matrix of $f(\cdot )$. Specifically, if all eigenvalues of the Jacobian evaluated at $x^{(0)}$, $\alpha ^{(0)}$ have negative real values, the steady state is asymptotically stable. If an eigenvalue crosses over into the right half of the complex plane as we move on the steady state manifold, the system becomes unstable. The critical boundary can therefore be characterized by the existence of an eigenvalue or a complex conjugate pair of eigenvalues on the imaginary axis. We state a system of equations $G(\cdot )= 0$ that formalizes this idea in (9) below.

Figure 4 shows the projection of the critical manifold on the parameter space. Stability can be ensured by enforcing parameter values on the stable side of the critical manifold. Robustness can be ensured by staying at a safe distance d from the critical manifold, where d must be chosen to account for the parametric uncertainties. More precisely, let $\alpha _i\in [\alpha _{i}^{(0)}-\Delta \alpha _i, \alpha _i^{(0)}+\Delta \alpha _i]$ represent any of the uncertain parameters stated in Table 2. After rescaling according to ${\tilde{\alpha }}_i= \alpha _i/\Delta \alpha _i$ the intervals

$$\begin{aligned} {\tilde{\alpha }}_i\in [{\tilde{\alpha }}_{i}^{(0)}-1, {\tilde{\alpha }}_i^{(0)}+1],\, i= 1, \dots , n_\alpha \end{aligned}$$

(7)

define a hypersquare, which is sketched for $n_\alpha = 2$ uncertain parameters in Fig. 4. Robustness is guaranteed by enforcing the center of this hypersquare to be at a distance larger than $\sqrt{n_\alpha }$ from the critical boundary (see Fig. 4) along the normal to the critical boundary. Formally, this results in the constraint

$$\begin{aligned} d-\sqrt{n_\alpha }\ge 0. \end{aligned}$$

The normal vector to the critical manifold is denoted by r.

Fig. 4

The normal vector r connects the parameter vector $\alpha ^{(0)}$ and the closest point $\alpha ^{(c)}$ on the critical manifold defined by $G(\cdot )= 0$. The square represents the uncertainty hypersquare (7). We assume $\Delta \alpha _i= 1$, or equivalently, ${\tilde{\alpha }}_i= \alpha _i$ for all i in the figure

It remains to state the equations $G(\cdot )= 0$. We need to characterize critical manifolds due to fold bifurcations, i.e., single real eigenvalues crossing the imaginary axis, because instability results from this case in the fluid catalytic cracker treated here. More generally, we treat critical manifolds that are characterized by the existence of an eigenvalue with real part $\sigma \le 0$. The case $\sigma = 0$ corresponds to the fold bifurcation. The case $\sigma < 0$ is used to enforce exponential stability with a prescribed decay rate $\sigma$ to ensure a sufficiently fast disturbance rejection. We call the critical manifold that corresponds to the latter case “modified fold manifold” for brevity.

The equations $G({\tilde{x}}^{(c)},\alpha ^{(c)},p, w)=0$, which have been abbreviated by $G(\cdot )= 0$ so far, read

$$\begin{aligned} f({\tilde{x}}^{(c)},{\tilde{x}}^{(c)},\ldots ,{\tilde{x}}^{(c)},\alpha ^{(c)},p)&=0 \end{aligned}$$

(8a)

$$\begin{aligned} \sigma w-A_0\, w-\sum _{i=1}^m A_i \exp (-\sigma \tau _i)\,w&=0 \end{aligned}$$

(8b)

$$\begin{aligned} w^\prime w\,-\,1&=0 . \end{aligned}$$

(8c)

where (8a) enforces ${\tilde{x}}^{(c)}$ to be a steady state for the parameters $\alpha ^{(c)}$. Equations (8b) ensure an eigenvalue $\lambda =\sigma$ to exist, where w is the corresponding eigenvector and the matrices $A_0$ and $A_i$ are the Jacobians of $f(x(t),x(t-\tau _1),...,x(t-\tau _m),\alpha ^{(c)},p)$ with respect to x(t) and $x(t-\tau _i)$, respectively, evaluated at $x(t)=x(t-\tau _i)={\tilde{x}}^{(c)}$. Because (8b) defines the direction of w but not its length, equation (8c) is required.

The system of equations that defines the normal vector r can be derived from (8) with a procedure first described in Mönnigmann and Marquardt (2002). We omit details of the derivation and only state the resulting system. It reads

$$\begin{aligned} G({\tilde{x}}^{(c)},\alpha ^{(c)},p, w)=0 \end{aligned}$$

(9a)

$$\begin{aligned} \begin{bmatrix} \nabla _{{\tilde{x}}^{(c)}}f^\prime &{}B_{12}&{}0\\ 0&{}B_{22}&{}2w&{} 0\\ \end{bmatrix}\kappa&= 0 \end{aligned}$$

(9b)

$$\begin{aligned} \begin{bmatrix} \nabla _{\alpha ^{(c)}}f^\prime &{} B_{32}&{}0\\ \end{bmatrix}\kappa -r&= 0 \end{aligned}$$

(9c)

$$\begin{aligned} r^\prime r- 1&=0 \, , \end{aligned}$$

(9d)

where

$$\begin{aligned} \begin{aligned} B^{{\mathrm {fold}}}_{12}&=-\sum _{i=0}^m\exp (-\sigma \tau _i)\nabla _{{\tilde{x}}^{(c)}}(w^\prime A_i^\prime )\\&\quad +\sigma \sum _{i=0}^m w^\prime A_i ^\prime \exp (-\sigma \tau _i)\nabla _{{\tilde{x}}^{(c)}}\tau _i\,,\\ B^{{\mathrm {fold}}}_{22}&= \sigma I-\sum _{i=0}^mA_i^\prime \exp (-\sigma \tau _i)\,,\\ B^{{\mathrm {fold}}}_{32}&=-\sum _{i=0}^m\exp (-\sigma \tau _i)\nabla _{\alpha ^\text {(c)}}(w^\prime A_i^\prime )\\&\quad +\sigma \sum _{i=0}^m (\nabla _{\alpha ^{\text {(c)}}}\tau _i)\big (\exp (-\sigma \tau _i)\,w^\prime \big )A_i^\prime \end{aligned} \end{aligned}$$

(10)

are the transposed Jacobians of (8b) w.r.t. ${\tilde{x}}^{(c)}$, w and $\alpha ^{(c)}$, respectively. The expression $\nabla _{ {\tilde{x}}^{(c)}} (w^\prime A_i^\prime )$ is given by

$$\begin{aligned} (\nabla _{ {\tilde{x}}^{(c)}} (w^\prime A_i^\prime ))_{\mu ,\nu }=\sum _{\rho =1}^n w_\rho \frac{\partial ^2 f_\nu }{\partial x _\mu \,\partial x_\rho (t-\tau _i)} \end{aligned}$$

(11)

and the expression $\nabla _{ \alpha ^{(c)}} (w^\prime A_i^\prime )$ is defined accordingly. Essentially, (9b) and (9c) span the normal space to the critical manifold in the combined space of the state variables and parameters and select the particular normal vector that has only components in the parameter space. Equation (9d) is required to normalize r to unit length (Mönnigmann and Marquardt 2002; Otten 2020). We abbreviate (9b)–(9d) by $H({\tilde{x}}^{(c)}, \alpha ^{(c)}, \kappa , r)= 0$ or $H(\cdot )= 0$. All derivatives in (9), (10) and (11) are calculated symbolically (Otten and Mönnigmann 2018).

Having defined the normal vector r, the robustness constraints can now be stated as

$$\begin{aligned} \alpha ^{(0)}-\alpha ^{(c)} + d\frac{r}{||r||}&= 0 \end{aligned}$$

(12a)

$$\begin{aligned} d- \sqrt{n_\alpha }&\ge 0. \end{aligned}$$

(12b)

We briefly note that optimization with the normal vector constraints need to be initialized at a, or close to a, feasible point. An algorithm for finding such initial points is described in (Otten and Mönnigmann 2016).

5 Optimization with robust stability constraints

Extending the optimization problem (6) by the normal vector constraints for robust stability results in

$$\begin{aligned}&\max \phi \end{aligned}$$

(13a)

$$\begin{aligned}&\quad \text {s. t.}\, f(x^{(0)},x^{(0)},x^{(0)},\alpha ^{(0)},p) = 0 \end{aligned}$$

(13b)

$$\begin{aligned}&\alpha ^{(0)}\text {within feasible intervals} \end{aligned}$$

(13c)

$$\begin{aligned}&2\,{\mathrm {kg}}/{\mathrm {s}}\le F_s\le 200\,{\mathrm {kg}}/{\mathrm {s}} \end{aligned}$$

(13d)

$$\begin{aligned}&300\,{\mathrm {K}}\le T_\text {air}\le 2500\,{\mathrm {K}} \end{aligned}$$

(13e)

$$\begin{aligned}&G(x^{(c)},\alpha ^{(c)},p,w) = 0 \end{aligned}$$

(13f)

$$\begin{aligned}&H(x^{(c)},\alpha ^{(c)},p,w,\kappa ,r) = 0 \end{aligned}$$

(13g)

$$\begin{aligned}&\alpha ^{(0)}-\alpha ^{(c)} + d\frac{r}{||r||} = 0 \end{aligned}$$

(13h)

$$\begin{aligned}&d\,-\sqrt{n_\alpha } \ge 0\, . \end{aligned}$$

(13i)

where (13a) to (13e) are equal to (6). The optimization is carried out with respect to $x^{(0)}$, $\alpha ^{(0)}$, p, $x^{(c)}$, $\alpha ^{(c)}$, w, r, $\kappa$ and d, where $F_s$ and $T_\text {air}$ can be removed by substituting (3). The optimal state variables and parameters are denoted by $x^{(0)}$ and $\alpha ^{(0)}$ to distinguish them from the closest critical point $x^{(c)}$, $\alpha ^{(c)}$. The remaining constraints comprise the augmented system $G(\cdot )= 0$ from (8) for the closest critical point $x^{(c)}$, $\alpha ^{(c)}$ in (13f), the normal vector system $H(\cdot )= 0$ from (9) that defines r at this critical point in (13g), and the distance constraints explained in (12). Optimization problems with normal vector constraints like (13) have been solved for systems modeled with ordinary differential equations before (Mönnigmann and Marquardt 2002; Mönnigmann et al. 2007). The delay differential equations (4) require a different augmented system (8) and normal vector system (9) that take the exponential function in the characteristic equation (8b) into account.

We compare results obtained without stability constraints and results obtained with the two stability requirements introduced in Sect. 4. The latter two are obtained by solving (13) for two different choices of $G(\cdot )= 0$, $H(\cdot )= 0$ in (13f), (13g). Specifically, we first set (13f), (13g) to the systems for fold bifurcation manifolds to ensure robust asymptotic stability. In the second case, we set (13f), (13g) to the systems for modified fold manifolds to guarantee robust exponential stability with a prescribed convergence rate to the steady state. We choose the convergence rate to correspond to a time constant of 10 min. This corresponds to

$$\begin{aligned} \sigma =-1/600\,{\mathrm {s}}^{-1}. \end{aligned}$$

(14)

The case without stability constraints corresponds to solving (6). It results in an unstable point of operation and is treated here for comparison only. The results for (6) and the two cases of (13) are listed in Table 4.

Before discussing the optimization results, problem (13) deserves some comments. We stress that problem (13) guarantees robust dynamic properties, such as asymptotic stability or a prescribed disturbance rejection rate in spite of uncertain parameters, without involving differential equations in the problem statement and solution procedure. By virtue of the normal vector method, the desired dynamic properties are achieved by enforcing a parametric distance of the optimal point to the closest critical boundary (cf. Fig. 4 and Section 4). This conversion of the dynamic problem into an algebraic one is particularly attractive for stiff differential equations, because a special numerical treatment of separated time scales is no longer required. The price to be paid for converting a dynamic problem into an algebraic one is the size of the resulting optimization problem (13). Specifically, (13b)–(13e) repeat the 95 equality constraints (6b), 14 inequality constraints (6c) and 4 inequality constraints (6d)–(6e) that already appeared in (6). The constraints (13f)–(13i) are 392 equality constraints and 1 inequality constraint, respectively, that constitute the normal vector constraints for robustness explained in Sect. 4. The number of optimization variables increases from 102 in (6) (95 state variables and the 7 parameters from Table 2) to 496 (additional variables $x^{(c)}$, $\alpha ^{(c)}$, w, r, $\kappa$ and d). In spite of the increased size, however, problem (13) can easily be solved with standard algorithms and implementations.¹

Table 4 FCC optimization results (see Tables 1 and 2 for units). The column labeled ’without’ corresponds to the optimization without normal vector constraints (6)

When comparing the optimization results, a profit reduction caused by stability requirements can be interpreted as the cost of stability. While an unstable point of operation results for (6) without the normal vector constraints, the profit is practically identical for this unstable optimum and the robust optima that result from (13) with asymptotic and exponential stability requirements (cf. Table 4). The cost of asymptotic as well as exponential stability with a prescribed decay rate of $\sigma =-1/600\,{\mathrm {s}}^{-1}<0$ is therefore negligible.

The optimization yields optimal feed temperatures and mass flow rates that are very similar in all three cases, while the parameters defining the temperature control for the regenerator and riser differ. The normal vector constraint (13i) enforcing robust asymptotic and robust exponential stability with the prescribed decay rate is active in both cases with normal vector constraints. This is consistent with the optimization without stability constraints resulting in an unstable steady state. Furthermore, the upper bounds on the catalyst circulation rate ($F_s\le 200\,{\mathrm {kg}}/{\mathrm {s}}$) and the gas oil feed ($F_{gR}\le 30\,{\mathrm {kg}}/{\mathrm {s}}$) are active for all three optimal points, the unstable one that results from (6) and those from (13) with constraints for asymptotic and exponential stability. It is plausible that these two constraints are active. As the gas oil feed $F_{gR}$ determines the amount of available reactant, it immediately affects the productivity of the FCC. The heat energy from the regeneration reactions can be transferred to the riser only via hot catalyst. The catalyst circulation rate $F_s$ therefore determines how much heat energy is available to power the endothermic cracking reactions.

The critical manifolds in the vicinity of the optimal points are shown in Figs. 5 and 6. These plots have to be interpreted with care, since they involve projections and cuts of a higher-dimensional parameter space. Figure 5 shows the optima in the $T_{\text {Reg,SP}}$-$F_\text {air}$-plane. Both the critical boundary for asymptotic stability due to fold bifurcations and the critical boundary for exponential stability due to modified fold points are shown. The corresponding two optimal points are located close to each other, which is also obvious from the numerical values given in Table 4. The hyperspherical outer approximations of the uncertainty regions touch the respecitve critical manifold, which confirms the constraints for robust stability are active (i.e., inequality (13i) is active). It is evident that the region of parametric uncertainty is located on the desired side of the critical manifold in both cases.

Fig. 5

Optimal points in $T_{\text {Reg,SP}}$-$F_\text {air}$-plane. Optimum and uncertainty region for asymptotic (dashed) and exponential stability (solid). Contour lines are shown for parameters set to their values at the optimum with robust asymptotic stability. a Overview of T_Reg,SP-Fair-plane with optima and critical manifolds. The region marked by the dotted rectangle is enlarged in (b). b Enlargement of the region marked by the dotted rectangle in (a)

Figure 6 shows the optima in the $T_{gR}$-$F_{gR}$-plane with the active upper bound on $F_{gR}$ and the critical manifolds. We stress that the critical manifolds and the hypercubical uncertainty regions only apparently overlap, because the critical manifolds are shown as cuts located at the critical point $\alpha ^{(c)}$, while the remaining features are shown as projections to the $T_{gR}$-$F_{gR}$-plane. Moreover, we stress that the upper boundary on $F_{gR}$ (horizontal line with grey hatched area) and the lower and upper boundaries on $T_{gR}$ (vertical lines with grey hatched areas) are shown such that a touching hypersphere corresponds to an active constraint. The upper bound on $F_{gR}$, for example, reads $F_{gR,\text {max}}= 30$kg/s according to Table 2. The corresponding constraint reads

$$\begin{aligned} F_{gR}\le F_{gR,\text {max}}+ \sqrt{n_\alpha }\Delta F_{gR}= 30\text {kg/s}+ \sqrt{5}\cdot 2\text {kg/s}= 34.47\text {kg/s} \end{aligned}$$

(15)

after overestimation with the Coke burning reaction rate coeicientf$(n_\alpha = 5)$-dimensional uncertainty hypersphere introduced in Fig. 4 and (12). The bound $F_{gR}\le 34.47$kg/s is visualized in Fig. 6. We choose to visualize (15) in Figs. 5 and 6, because an active constraint corresponds to a tangential intersection of the uncertainty hypersphere and a critical boundary for all boundaries in this case. If the actual constraint was visualized instead of (15, a touching hypersphere would correspond to an active constraint for stability boundaries, while a touching of its center point would correspond to an active feasibility constraint. We stress again this is only a matter of visualization and the original bounds and the ones modified as in (15) are equivalent. The corresponding boundaries for $T_{gR}$ in Fig. 5 are obtained accordingly, which results in $T_{gR}\le 600\text {K}+ \sqrt{5}\cdot 20\text {K}= 644.7\text {K}$ and $T_{gR}\ge 300\text {K}- \sqrt{5}\cdot 20\text {K}= 255.3\text {K}$.

Fig. 6

Optimal points in the $T_{gR}$-$F_{gR}$-plane. Optimum and uncertainty region for asymptotic stability (dashed) and exponential stability (solid). Contour lines are shown for parameters set to their values at the optimum with robust asymptotic stability. a Overview of T_gR-F_gR-plane with optima and critical manifolds. The region marked by the dotted rectangle is enlarged in (b). b Enlargement of the region marked by the dotted rectangle in (a)

We corroborate the robustness of the optima that result from (13) with simulations. We select the optimum obtained with the normal vector constraints for exponential stability (right column in Table 4) for these illustrations. The corresponding results for the other robust optimum (center column in Table 4) are very similar and omitted for brevity. We pointed out earlier that the system model consists of stiff delay differential equations, because the riser dynamics are much faster than the regenerator dynamics. We use a toolbox tailored to the numerical integration of stiff delay differential equations for all simulations (Rackauckas and Nie 2017).

Figure 7 shows a simulation of riser and regenerator temperatures after a step disturbance of the steady state. It is evident that the system is stable.

Fig. 7

Temperatures in the simulated FCC. The regenerator temperature is lowered by $10\%$ to disturb the FCC at $t=0$

A time series of concentrations in the regenerator gas phase is shown in Fig. 8. The reduced regenerator temperature initially inhibits the burning of coke, which results in an immediate rise in $y_{0_2}$ and an immediate drop in $y_{CO_2}$. The concentrations then approach their new steady state values.

Fig. 8

Simulation of gas phase in the regenerator after the same disturbance as in Fig. 7. The reduced temperature lowers the reaction rates, which causes the oxygen concentration (reactant) to rise and the carbon dioxide concentration (reaction product) to drop temporarily

Figure 9 shows the first five seconds of the same simulation. The initial disturbance has an impact on the riser temperature with a time delay because the cool catalyst from the regenerator has to pass a standpipe connecting regenerator and riser.

Fig. 9

Temperatures in the simulated FCC. Only a short time span after the simulations start is shown. The delayed riser temperature drop is caused by the transport delay due to standpipes connecting regenerator and riser

In Fig. 10, the concentrations $y_i$ are shown against the riser length at a steady state. It visualizes the intended productive reaction. The gas oil concentration $y_A$ drops as the gas oil is converted to gasoline, coke and light hydrocarbons, causing the concentrations $y_B$, $y_C$ and $y_D$ to increase while traveling through the riser.

Fig. 10

Components $y_i$ in the riser at steady state. The concentration of gas oil $y_A$ drops as it is converted to gasoline, coke and light hydrocarbons as expected

Finally, Fig. 11 illustrates the prescribed robustness has been achieved. The figure shows the time series for ten representative states, including those from Figs. 7 and 8 , with lower and upper bounds that correspond to the decay rate (14). It is evident that the bounds are respected.

Fig. 11

Time series of ten representative states that result after the regenerator temperature is artificially lowered by 10 % to cause a disturbance at $t = 0$. Solid lines mark trajectories, dashed lines mark the lower and upper bounds that result for the decay rate (14). $T_{\text {Reg}}$, $y_{O_2}$, $y_{CO}$ and $y_{CO_2}$ are regenerator states. The remaining variables belong to segment $l=15$ at the upper end of the spatially discretized riser

6 Conclusion

We showed the normal vector method for robust optimization of parametrically uncertain systems can be applied to stiff delay differential equations. The viability of the method was demonstrated by applying it to an FCC model that consists of 95 differential equations with two delays and five uncertain parameters. The method was successfully used to ensure robust asymptotic stability and robust exponential stability with a prescribed disturbance rejection rate.

Model parameters and constants

See Tables 5, 6 and 7.

Table 5 FCC model abbreviations and substitutions

Table 6 FCC model constants: Roman letters

Table 7 FCC model constants: greek etters