Handling Measurement Delay in Iterative Real-Time Optimization Methods

Abstract

This paper is concerned with the real-time optimization (RTO) of chemical plants, i.e., the optimization of the steady-state operating points during operation, based on inaccurate models. Specifically, modifier adaptation is employed to cope with the plant-model mismatch, which corrects the plant model and the constraint functions by bias and gradient correction terms that are computed from measured variables at the steady-states of the plant. This implies that the sampling time of the iterative RTO scheme is lower-bounded by the time to reach a new steady-state after the previously computed inputs were applied. If analytical process measurements (PAT technology) are used to obtain the steady-state responses, time delays occur due to the measurement delay of the PAT device and due to the transportation delay if the samples are transported to the instrument via pipes. This situation is quite common because the PAT devices can often only be installed at a certain distance from the measurement location. The presence of these time delays slows down the iterative real-time optimization, as the time from the application of a new set of inputs to receiving the steady-state information increases further. In this paper, a proactive perturbation scheme is proposed to efficiently utilize the idle time by intelligently scheduling the process inputs taking into account the time delays to obtain the steady-state process measurements. The performance of the proposed proactive perturbation scheme is demonstrated for two examples, the Williams–Otto reactor benchmark and a lithiation process. The simulation results show that the proposed proactive perturbation scheme can speed up the convergence to the true plant optimum significantly.

1. Introduction

There is a strong interest in the process industries to identify the best steady-state operating conditions of their processes such that they are operated at their economic optimum while meeting constraints on product qualities, emissions, and equipment capabilities. Optimal operating conditions of a process can be identified by formulating and solving an optimization problem where a cost (or profit) function is minimized (or maximized) respecting all process limitations on the basis of a mathematical model of the plant. As the parameters of the optimization problem (e.g., prices for raw materials, energy, and products) and the behavior of the plant vary over time, such optimization is performed repeatedly during plant operation, which is called RTO [1,2] (real-time optimization). At each call of the optimizer, the problem at hand is a constrained nonlinear program (NLP) which usually is solved by gradient-based algorithms, e.g., SQP solvers [3,4] or interior point algorithms [5,6]. These solvers compute local optima, which is sufficient under the tacit assumption that only minor adaptations are made and the problem is convex in the region of interest. As far as the model-based optimization is concerned, in principle also other optimization methods can be applied, as e.g., derivative free methods [7,8], population-based methods [9], or nature-inspired algorithms [10]. While for the solution of the optimization problems efficient algorithms are available, a major practical problem is the mismatch between the model and the true behavior of the plant. No matter how the optimum is computed, if the model and the behavior of the plant do not match, the computed optimum won’t be the true optimum for the plant and the solution may even violate the constraints of the plant. Of course the model could be discarded completely and a brute force optimization could be performed using any derivative-free or meta-heuristic algorithm, but this implies performing a large number of exploratory moves with the real plant, which is highly undesirable. In particular, large constraint violations of internal variables must be avoided, which is impossible if the optimization is completely model-free. Therefore, a combination of model-based and data-based optimization is advantageous. The classical approach to this problem is to add a parameter estimation element, called the two-stage approach [11].

In the two-step approach, some model parameters are updated repeatedly using least-squares estimation and the optimization problem with updated model parameters is solved to identify the process optimum. The performance of the two-step approach is discussed in detail in Chen and Joseph [12]. The two-step approach is designed to handle parametric plant-model mismatch, which means that the model can, with the correct parameter values, represent the true plant exactly. However, there always is also structural plant-model mismatch, i.e., even with optimized parameters, the predictions of the model are inaccurate. This issue was addressed in the “Integrated system optimization and parameter estimation” (ISOPE) approach proposed in Roberts [13]. In ISOPE scheme, in addition to updating the model parameters in each iteration, the objective function of the model based optimization problem is modified iteratively by adding bias and gradient correction terms so that the solution converges iteratively to the optimum of the real plant. These correction terms (also called modifiers) are estimated from the available measurements.

Later, by Tatjewski [14], it was shown that the model parameters do not have to be updated in each iteration to converge to a process optimum. Based on this insight, Tatjewski [14] simplified the ISOPE approach and proposed the redesigned-ISOPE scheme without the parameter-update step from the ISOPE approach. The redesigned-ISOPE scheme was extended in Gao and Engell [15] to handle process-dependent constraints, by adding bias and gradient correction terms also to the constraint functions, and was termed iterative gradient modification optimization (IGMO). The IGMO scheme in Gao and Engell [15] was analyzed in detail in Marchetti et al. [16] and there was given the name modifier adaptation (MA).

The efficiency of MA-based iterative RTO methods depends on the accuracy of the gradient modifiers, which contain the gradients of the response of the cost function and of the constraints to the inputs. As the plant gradients usually cannot be measured directly, they have to be computed from the available measurements. The traditional way of doing this is to use finite differences, preferably employing the past input moves as in Gao and Engell [15] where this was enhanced with probing moves if the gradient estimation problem becomes ill-conditioned. However, in finite differences-based schemes, a trade-off between the approximation error and the effect of measurement noise has to be found. Therefore, Gao et al. [17] proposed to use quadratic approximation (QA) for the approximation of the plant gradients, which is more robust at measuring noise compared to that of finite differences (FD) and Broyden’s formula [18]. The resulting modifier adaptation with quadratic approximation (MAWQA) scheme combines elements of MA, QA, and derivative free optimization [7] (DFO).

For all adaptation-based iterative RTO methods, it is a prerequisite that the adapted process model is adequate at the—unknown—optimum of the real plant [1,16,19,20,21]. A model is considered adequate if the Lagrangian of the model-based optimization problem satisfies the necessary second-order conditions of optimality at the optimum of the real process. In Faulwasser and Bonvin [22], it was proposed to use second-order correction terms (modifiers) in the objective and constraint functions to match the Hessian of the nominal model to that of the plant. However, it is difficult to estimate the plant Hessian [23] from real data. Model adequacy was addressed in François and Bonvin [24] by using convex model approximation of the objective and constraint functions. This may however slow down the rate of convergence to the process optimum. Ahmad et al. [25] addressed model adequacy by updating some model parameters only if the model adequacy conditions are satisfied at the current operating input. In Gottu Mukkula and Engell [26], a guaranteed model adequacy (GMA) scheme for MAWQA was proposed where model adequacy is ensured by strictly convex quadratic approximations of the objective and constraint functions if the goal of the optimization is to minimize an objective function. The MAWQA with GMA scheme may converge to the process optimum in fewer iterations than MAWQA [26].

No matter how the optimization model is corrected, by parameter adaptation or by modifier adaptation of by a combination of both [23], the solution of the resulting model-based optimization problem can be performed by any algorithm. In our work, we use gradient-based algorithms, as the cost functions and the constraints are smooth, no discrete variables are involved, and guaranteed (local) optimality is important.

Modifier adaptation-based iterative RTO methods require steady-state process measurements to compute a new input that provides a lower limit to the time between iterates. This can be alleviated by the use of transient process measurements [27,28,29,30,31], but this is outside the scope of this paper. In addition to the time required to reach a steady-state, the convergence to a process optimum may also be slowed down due to delays in obtaining measurement information. Considerable measurement delays may occur due to the time required for the measurement device (e.g., gas chromatography) to analyze the sample and due to the remote positioning of a measurement device. Figure 1 illustrates a situation when the time delay is caused by the remote positioning of the measurement device for safety reasons, e.g., due to harsh conditions at the plant. Measurement devices like NMR have to be placed in a certified (e.g., ATEX, IECEx) enclosure located at a distance from the process equipment. The thin, long tubings that carry the sample from plant to the remotely positioned measurement device can cause a significant delay, and this is quite common in the process industries.

Gottu Mukkula et al. [32,33] proposed two schemes where additional plant perturbations are performed during the waiting period. The purpose was to gain additional information about the process instead of remaining idle until the effect of a new operating input propagated through the measurement device setup. Here, we propose a strategy in the context of MA that proactively schedules the input moves and thereby significantly reduces the effect of the time delay caused due to the remote positioning of a measurement device. The input perturbations in the proposed method are computed by solving an optimization problem using the latest measurement information.

The rest of the paper is organized as follows. Firstly, the MA scheme used in this paper is presented in detail. Then, the previously proposed approaches to handle the time delay due to the positioning of a measurement device are presented. Thereafter, the new proactive perturbation scheme to effectively handle the time delay is described. Finally, the performance of the proposed scheme is analyzed using the Williams–Otto reactor benchmark and a lithiation reaction case study.

2. Preliminaries

2.1. Model

Consider the steady-state mathematical models $\mathit{y}={\mathit{F}}_p\left(\mathit{u}\right)$ and $\widehat{\mathit{y}}={\mathit{F}}_m\left(\mathit{u}\right)$ as the exact representation of a process and as the nominal model that was built to represent the underlying process. The mapping functions ${\mathit{F}}_p:{\mathbb{R}}^{n_u}\to {\mathbb{R}}^{n_y}$ and ${\mathit{F}}_m:{\mathbb{R}}^{n_u}\to {\mathbb{R}}^{n_y}$ map the $n_u$-dimensional vector of manipulated variables $\mathit{u}$ to the $n_y$-dimensional vector of measured variables $\mathit{y}$ and $\widehat{\mathit{y}}$. Let $J:{\mathbb{R}}^{n_u}\times {\mathbb{R}}^{n_y}\to \mathbb{R}$ represent the objective function, which is continuous and twice differentiable with respect to $\mathit{u}$, that should be minimized, and let $\mathit{G}:{\mathbb{R}}^{n_u}\times {\mathbb{R}}^{n_y}\to {\mathbb{R}}^{n_c}$ be a continuous and twice differentiable with respect to $\mathit{u}$ $n_c$-dimensional vector of constraint functions, which may include process, safety, and quality restrictions.

An optimal input of the process lying within ${\mathit{u}}^L$ and ${\mathit{u}}^U$ can be obtained theoretically by solving

$$\begin{array}{cc}\hfill {\mathit{u}}_p^{*}=arg& \underset{\mathit{u}\in [{\mathit{u}}^L,{\mathit{u}}^U]}{min}J(\mathit{y},\mathit{u})\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.\mathit{y}={\mathit{F}}_p\left(\mathit{u}\right),\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill & \mathit{G}(\mathit{y},\mathit{u})\le \mathbf{0},\hfill \end{array}$$

(1c)

but the mapping function (1b) is not available. The optimum ${\mathit{u}}_m^{*}$ results from solving the optimization problem using the nominal model:

$$\begin{array}{cc}\hfill {\mathit{u}}_m^{*}=arg& \underset{\mathit{u}\in [{\mathit{u}}^L,{\mathit{u}}^U]}{min}J(\widehat{\mathit{y}},\mathit{u})\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.\widehat{\mathit{y}}={\mathit{F}}_m\left(\mathit{u}\right),\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill & \mathit{G}(\widehat{\mathit{y}},\mathit{u})\le \mathbf{0},\hfill \end{array}$$

(2c)

and may differ from ${\mathit{u}}_p^{*}$ considerably if ${\mathit{F}}_m\ne {\mathit{F}}_p$. To simplify the notation, we from here on replace $J(\mathit{y},\mathit{u}),\mathit{G}(\mathit{y},\mathit{u})$ by $J_p\left(\mathit{u}\right),{\mathit{G}}_p\left(\mathit{u}\right)$ and $J(\widehat{\mathit{y}},\mathit{u}),\mathit{G}(\widehat{\mathit{y}},\mathit{u})$ by $J_m\left(\mathit{u}\right),{\mathit{G}}_m\left(\mathit{u}\right)$.

2.2. Modifier Adaptation

Modifier adaption handles plant-model mismatch by adding and iteratively updating gradient correction term to the objective function, and by adding and iteratively updating bias and gradient correction terms to the constraint functions of the model based optimization problem (2) as described in Gao and Engell [15], Marchetti et al. [16]. In the kth iteration, the modified objective and constraint functions of the model based optimization problem, i.e., the MA problem, are

$$\begin{array}{cc}\hfill J_m^{ad,k}\left(\mathit{u}\right)& :=J_m\left(\mathit{u}\right)+{(\nabla J_p^k-\nabla J_m^k)}^T(\mathit{u}-{\mathit{u}}^k),\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill {\mathit{G}}_m^{ad,k}\left(\mathit{u}\right)& :={\mathit{G}}_m\left(\mathit{u}\right)+({\mathit{G}}_p^k-{\mathit{G}}_m^k)+{(\nabla {\mathit{G}}_p^k-\nabla {\mathit{G}}_m^k)}^T(\mathit{u}-{\mathit{u}}^k),\hfill \end{array}$$

(3b)

where $\nabla J_p^k,\nabla J_m^k$ are the gradients of the plant objective function $J_p\left(\mathit{u}\right)$ and of the nominal model objective function $J_m\left(\mathit{u}\right)$ with respect to the process input for the kth iteration. Similarly, $\nabla {\mathit{G}}_p^k$, $\nabla {\mathit{G}}_m^k$ represent the gradients of the plant constraint functions ${\mathit{G}}_p\left(\mathit{u}\right)$ and the constraint functions of the nominal model ${\mathit{G}}_m\left(\mathit{u}\right)$ with respect to the process input for the kth iteration. The modified optimization problem that is solved in the kth iteration is:

$$\begin{array}{cc}\hfill {\widehat{\mathit{u}}}^{k+1}=arg& \underset{\mathit{u}\in [{\mathit{u}}^L,{\mathit{u}}^U]}{min}{J}_m^{ad,k}\left(\mathit{u}\right)\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.{\mathit{G}}_m^{ad,k}\left(\mathit{u}\right)\le \mathbf{0},\hfill \end{array}$$

(4b)

where the determined optimum $({\widehat{\mathit{u}}}^{k+1})$ is the $k+1$th input which is applied to the plant.

The gradients of the plant objective function ($\nabla J_p^k$) and of the plant constraint functions ($\nabla {\mathit{G}}_p^k$) can be obtained for example using finite differences for which additional plant perturbations are required. Methods to approximate the plant gradients using the available past ℓ number of inputs and steady-state measurements, without the need for additional plant perturbations, are proposed in Roberts [18], Brdyś and Tatjewski [34], and Gao et al. [35].

In MAWQA, the gradients of the plant objective function ($\nabla J_p^k$) and of the plant constraint functions ($\nabla {\mathit{G}}_p^k$) are obtained using surrogate QA models as long as the cardinality condition, i.e., the minimum number of data points required to identify all the parameters of the QA function, is fulfilled. A minimum of $\ell :=\frac{(n_u+1)(n_u+2)}{2}$ data points, i.e., process inputs and their corresponding steady-state process measurements, are required to fit a surrogate QA function (Q) for the plant objective function and for each of the plant constraint functions. Upon satisfying the cardinality condition, $\nabla J_p^k$ and $\nabla {\mathit{G}}_p^k$ are then obtained by evaluating the analytical derivatives of the surrogate QA models at ${\mathit{u}}^k$. A surrogate QA model can be defined as

$$Q(\mathit{p},\mathit{u})=\sum _{i=1}^{n_u}\sum _{j=1}^i{a}_{i,j}{u}_i{u}_j+\sum _{i=1}^{n_u}{b}_i{u}_i+c,$$

(5)

where, $\mathit{p}:=\{a_{1,1},\dots ,a_{n_u,n_u},b_1,\dots ,b_{n_u},c\}$ are the parameters of the surrogate QA model.

To obtain a surrogate QA model for each of the objective and constraint functions, a subset of at least ℓ data points (${\mathcal{U}}^k$) is selected from the set of all data points available until the kth iteration (${\mathbb{U}}^k$). Wenzel et al. [36] reported a comparative study of existing methods to screen ${\mathcal{U}}^k$ from ${\mathbb{U}}^k$. In general, the methods for screening ${\mathbb{U}}^k$ attempt to include distant data points which are well distributed around ${\mathit{u}}^k$ (${\mathcal{U}}_{dist}^k$) and can act as anchor points for the QA. They also include all neighboring data points (${\mathcal{U}}_{nb}^k$), which lie within the inner-circle around ${\mathit{u}}^k$. The inner-circle is defined as an n-dimensional sphere around ${\mathit{u}}^k$ with the tuning parameter Δu as its radius. As proposed in Gao et al. [17], the quality of the distribution of data points in ${\mathcal{U}}_{dist}^k$ can be evaluated using a so-called conditionality-check criterion, the inverse of the condition number of ${\mathit{s}}^k$ (${\kappa }^{-1}\left({\mathit{s}}^k\right)$) has to be greater than a desired value ($\delta$). ${\mathit{s}}^k$ in the conditionality-check criterion is defined as

$${\mathit{s}}^k={[{\mathit{u}}^k]}^{n_u\times 1}\otimes 1^{\mathrm{cardinality}\left({\mathcal{U}}_{dist}^k\right)\times 1}-[{\mathcal{U}}_{dist}^k],$$

(6)

where $[{\mathcal{U}}_{dist}^k]$ and $[{\mathit{u}}^k]$ are matrix representations of the set ${\mathcal{U}}_{dist}^k$ and of ${\mathit{u}}^k$. If the conditionality-check criterion fails, additional plant perturbations are performed and added to ${\mathcal{U}}^k$ to improve the distribution of the anchor points used for fitting the QA-function.

As the data points in ${\mathcal{U}}^k$ consist of input variables positioned locally around ${\mathit{u}}^k$, the QA functions fitted using ${\mathcal{U}}^k$ are also local approximations of $J_p$ and ${\mathit{G}}_p$. In MAWQA, this is taken into account by restricting the process input for the next iteration (${\widehat{\mathit{u}}}^{k+1}$) obtained by solving the optimization (4) to lie inside an ellipsoidal trust region by adding the following constraint to the modifier adaptation problem in (4)

$${(\mathit{u}-{\mathit{u}}^k)}^T\mathrm{cov}\left({\mathcal{U}}^k\right)(\mathit{u}-{\mathit{u}}^k)\le {\gamma }^2,$$

(7)

where $\gamma$ is a tuning parameter that scales the trust region [17]. In addition, there are also elements of DFO [7] that include the criticality-check, the quality-check and a possible switch to an optimization based on the quadratic approximation model. We refer the readers to Gao et al. [17,37] for a detailed discussion on MAWQA.

Convergence to the plant optimum is assumed when the evaluation function $\mathcal{C}$ defined as ${\mathcal{C}}_i\stackrel{\mathsf{def}}{=}{\sum }_{j=1}^{n_u}\left(\right|u_j^i-u_j^{i-1}\left|\right|\nabla J_{p,j}^i\left|\right)$ is less than $|\epsilon J_{p,best}^k|\forall i=\{k+1-N_c,\dots ,k\}$ where $\nabla J_{p,j}^i$ is the gradient of the plant objective function of the jth input-variable of kth input $\left({\mathit{u}}^k\right)$ [38]. The tuning parameter $N_c$ is the number of iterations for which the condition $\mathcal{C}\le |\epsilon J_{p,best}^k|$ has to be satisfied continuously. $J_{p,best}^k$ is the best value of the objective function realized by the plant up to the kth-iteration and $\epsilon$ is a tuning parameter that must be chosen taking into account the measurement noise and the required tolerance of the value of the optimum.

3. Active Perturbation Strategies

Figure 1 illustrates a situation that occurs typically in the process industries while implementing MA based RTO methods to optimize their processes. Consider a process that reaches a steady-state for constant inputs within ${\tau }_p$ time units and the measurement device that requires a maximum of ${\tau }_m$ time units to analyze a sample. An additional ${\tau }_d$ time units is required for the sample to be transported to the measurement device making the total time to obtain steady-state measurements after a change of the inputs ${\tau }_t:={\tau }_p+{\tau }_d+{\tau }_m$. If the measurement device is located at a remote position, this leads to a waiting period.

Figure 2 illustrates the waiting period in an iteration of a MA algorithm. If an input $u^{k+1}$ is given to the plant at time t, its corresponding steady-state plant measurement becomes available at time $t+{\tau }_t$. The time between iterations in MA is thus at least ${\tau }_t$ time units if the plant gradients are computed using past inputs and measurements, and even larger if additional perturbations are performed. During the waiting period the plant remains at a steady-state and no additional process information is gained during this period, hence the convergence to the optimum steady-state is slowed down. Gottu Mukkula et al. [32,33] proposed active perturbation strategies to perform input perturbations during the waiting period to gain additional process information that is used later to drive the process to its optimum faster. Figure 2 illustrates the active perturbation scheme where additional process perturbing inputs are given to the process, for example, between $t+{\tau }_t+{\tau }_p$ and $t+2{\tau }_t$, to gain additional process information that can be used in the subsequent iterations, once the measurements of the effects of the input perturbations arrive. In general, the maximum possible number of input perturbations that can be given to a process within the waiting period can be computed as [32,33]

$$P_{max}:=\lfloor min\left\{\frac{{\tau }_d+{\tau }_m}{{\tau }_p},\frac{{\tau }_d+{\tau }_p}{{\tau }_m}\right\}\rfloor .$$

(8)

In the following subsections, we present the active perturbation strategies proposed in Gottu Mukkula et al. [32,33]. Both active perturbation strategies are illustrated using an example case with 2 input variables and by assuming the maximum possible perturbations per iteration ($P_{max}$) as 4.

3.1. Active Perturbation around the Current Input (APCI)

According to the method proposed in Gottu Mukkula et al. [32], perturbing inputs around the input in the current iteration, i.e., ${\mathit{u}}^k$, are given to the process during the waiting period. Figure 3, illustrates this active perturbation method. In the figure, represents the inputs from MA. Let ${\mathit{u}}^k$ be the input given to the plant in the kth iteration at time t. The plant reaches its steady-state at $t+{\tau }_p$, and the steady-state measurements for ${\mathit{u}}^k$ are available at $t+{\tau }_t$. During the waiting period, i.e., from $t+{\tau }_p$ to $t+{\tau }_t$, $P_{max}=4$ additional perturbations represented by circles (◯) are performed around ${\mathit{u}}^k$. Redundant probing is avoided by suppressing perturbation inputs that are closer than a threshold $\xi$ from one of the earlier inputs. For example, the perturbing input represented by $\oslash$ is suppressed as it is in close vicinity to ${\mathit{u}}^{k-1}$, and therefore, would be redundant. During the remaining waiting period, which is available due to suppressing of some perturbations, the plant is given the best known input $({\mathit{u}}_{p,best}^k)$ to avoid a loss of performance. As the perturbing inputs are not computed optimally, not all of them provide useful information. For example, in Figure 3, only the perturbation inputs represented by $\oplus$ are used by the QA scheme for computing the gradient at ${\mathit{u}}^{k+1}$.

3.2. Active Perturbation around an Estimate of the Next Input (APENI)

Here, the process is perturbed around an estimate of the next input, i.e., ${\widehat{\mathit{u}}}_e^{k+1}$, instead of perturbing around current input ${\mathit{u}}^k$. Due to the unavailability of the measurements for the input ${\mathit{u}}^k$ during the waiting period, the plant gradients at ${\mathit{u}}^k$ cannot not be computed. So, the next input ${\mathit{u}}^{k+1}$ cannot be calculated precisely without the plant gradients at ${\mathit{u}}^k$, and therefore is estimated using the following assumption: If ${\mathit{u}}^k$ is close to ${\mathit{u}}^{k-1}$, the difference between the plant and the model gradients at ${\mathit{u}}^k$ is approximately equal to the gradient difference at ${\mathit{u}}^{k-1}$, i.e., $(\nabla J_p^k-\nabla J_m^k)\approx (\nabla J_p^{k-1}-\nabla J_m^{k-1})$ and $(\nabla {\mathit{G}}_p^k-\nabla {\mathit{G}}_m^k)\approx (\nabla {\mathit{G}}_p^{k-1}-\nabla {\mathit{G}}_m^{k-1})$. The same is assumed for the bias correction term $({\mathit{G}}_p^k-{\mathit{G}}_m^k)$. In other words, the optimization problem is unchanged. The input in the next iteration can then be estimated by solving the following modified MA problem:

$$\begin{array}{cc}\hfill {\widehat{\mathit{u}}}_e^{k+1}=arg& \underset{\mathit{u}\in [{\mathit{u}}^L,{\mathit{u}}^U]}{min}{J}_{m,e}^{ad,k}\left(\mathit{u}\right)\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.{\mathit{G}}_{m,e}^{ad,k}\left(\mathit{u}\right)\le \mathbf{0},\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill & {(\mathit{u}-{\mathit{u}}^k)}^T\mathrm{cov}\left({\mathcal{U}}^k\right)(\mathit{u}-{\mathit{u}}^k)\le {\gamma }^2,\hfill \end{array}$$

(9c)

where

$$\begin{array}{cc}\hfill J_{m,e}^{ad,k}\left(\mathit{u}\right)& =\begin{array}{cc}& J_m\left(\mathit{u}\right)+{(\nabla J_p^{k-1}-\nabla J_m^{k-1})}^T(\mathit{u}-{\mathit{u}}^k),\hfill \end{array}\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\mathit{G}}_{m,e}^{ad,k}\left(\mathit{u}\right)& =\begin{array}{cc}& {\mathit{G}}_m\left(\mathit{u}\right)+({\mathit{G}}_p^{k-1}-{\mathit{G}}_m^{k-1})+{(\nabla {\mathit{G}}_p^{k-1}-\nabla {\mathit{G}}_m^{k-1})}^T(\mathit{u}-{\mathit{u}}^{k-1}).\hfill \end{array}\hfill \end{array}$$

(10b)

The problem in (9) for computing ${\widehat{\mathit{u}}}_e^{k+1}$ is the same as the problem in (3) for computing ${\mathit{u}}^{k+1}$, except that the trust-region in (7) for computing ${\widehat{\mathit{u}}}_e^{k+1}$ is now centered around ${\mathit{u}}^k$.

In Figure 4a, an ideal scenario for an active perturbation around next input is illustrated. In the ideal case, the assumption on the modifiers is valid, i.e., $(\nabla J_p^k-\nabla J_m^k)\approx (\nabla J_p^{k-1}-\nabla J_m^{k-1})$ and $(\nabla {\mathit{G}}_p^k-\nabla {\mathit{G}}_m^k)\approx (\nabla {\mathit{G}}_p^{k-1}-\nabla {\mathit{G}}_m^{k-1})$. Then, ${\widehat{\mathit{u}}}_e^{k+1}\approx {\mathit{u}}^{k+1}$, and therefore, the measurement information from the perturbed inputs are used to compute the inputs in upcoming iterations. Figure 4b illustrates a nonideal scenario where the assumption about the modifiers does not hold, so ${\widehat{\mathit{u}}}_e^{k+1}\ne {\mathit{u}}^{k+1}$. Therefore, as shown in Figure 4b, not all measurement information from the perturbing inputs may be useful in the upcoming iterations. Redundant perturbation inputs are suppressed similar to the active perturbation around the current input (APCI) method.

4. Proactive Perturbation Scheme

Although the perturbation schemes proposed in Gottu Mukkula et al. [32,33] gain additional measurement information by perturbing the process during the waiting period, they are suboptimal as they are based on heuristics. Also, the measurement information gained from the additional input perturbations is not used immediately. Therefore, in this paper, we propose a proactive perturbation scheme where the additional perturbation inputs are computed by solving an optimization problem using all available process information. The key idea of the proactive perturbation scheme is illustrated in Figure 5 using an example case with $n_u=1$, $P_{max}=5$, and where the plant gradients are estimated using QA. At $t:=0$, an input $u^0$ is given to the process, the corresponding measurements become available after time ${\tau }_t$. In the initial (0th) iteration, FD is used for gradient approximation as there are no past data points for gradient approximation. Therefore, after applying $u^0$ for $max\{{\tau }_p,{\tau }_m\}$ units of time, the process is operated with $u_{f{d}_1}^0$, for gradient approximation, the corresponding measurements are available at time $\zeta :={\tau }_t+(n_u\times max\left\{{\tau }_p,{\tau }_m\right\})$ and a new iteration input $u^6$ is computed by solving the MA problem in (3). During the waiting period in the 0th iteration, i.e., from $(n_u+1)\times max\left\{{\tau }_p,{\tau }_m\right\}$ to $\zeta$, $P_{max}$ additional perturbation moves (from $u^1$ to $u^5$) around $u^0$ are performed to gain additional knowledge about the process. For $n_u:=1$, steady-state measurements for at least $\ell :=3$ different inputs are required for the process gradients to be approximated using QA (cardinality condition). We note that ℓ varies if other methods are used for plant gradient estimation using past data. In this example, after time $\zeta +max\left\{{\tau }_p,{\tau }_m\right\}$, there are already measurements for ℓ inputs. From here on, an MA problem is solved each time a new measurement is available as the plant gradients are computed using past inputs and measurements. From this point onward, each input to the process is optimal with respect to the available process knowledge, thereby making the proposed proactive perturbation scheme more efficient than the active perturbation methods proposed in Gottu Mukkula et al. [32,33].

A flowsheet for the implementation of the proposed proactive perturbation scheme is presented in Figure 6. In the flowsheet, $\varphi$ is the total number of inputs given to the plant, also including the inputs whose steady-state measurements are not yet available and $\psi$ is the total number of inputs given to the plant for whom the steady-state plant measurements are already available. The major steps in the flowsheet are:

➀

perturbation input for gradient estimation using FD

➁

perturbation input for gradient estimation using FD, also considered as an input (${\mathit{u}}^k$) for MA problem formulation

➂

perturbation input to gain additional information during the waiting period

➃

perturbation input to gain additional information during the waiting period, also considered as an input (${\mathit{u}}^k$) for MA problem formulation

➄

input (${\mathit{u}}^{k+1}$) computed by solving a MA problem (4) with plant gradients approximated using FD

➅

perturbation input considered as an input (${\mathit{u}}^k$) for MA problem formulation, to ensure a continuous flow of steady-state measurements every $max\left\{{\tau }_p,{\tau }_m\right\}$ units of time

➆

input (${\mathit{u}}^{k+1}$) computed by solving MA problem with gradients approximated from past (at least ℓ) measurements

In all the active perturbation strategies discussed in this section, to avoid frequent probing of the plant with unsuccessful inputs, and thus, deteriorating the plant performance, a list of unsuccessful moves is stored. An input from MA ${\mathit{u}}^k$ is applied to the plant only if the number of past inputs on the list that are not farther than $\phi$ from ${\mathit{u}}^k$ is less than a predefined number ${\chi }_{count}$. The parameters $\phi ,{\chi }_{count}$ are defined depending on the expected level of measurement noise.

5. Williams–Otto Reactor Case Study

5.1. Process Description

The Williams–Otto reactor case study [39] is a benchmark case study widely used to analyze the performance of RTO methods. Reactants A and B react in a CSTR to produce products E and P, and a side product G. While the reaction in the plant follows a three-step reaction mechanism, the nominal model considers a two-step reaction mechanism, ignoring the formation of an intermediate component C, leading to a structural plant-model mismatch. The reaction mechanisms for the plant and in the nominal model are:

$$\begin{array}{cccc}\hfill & \mathrm{Plant}:\hfill & \hfill & \mathrm{Nominal}\mathrm{model}:\hfill \\ \hfill & \mathrm{A}+\mathrm{B}\stackrel{k_1}{⟶}\mathrm{C},\hfill & \hfill & \mathrm{A}+2\mathrm{B}\stackrel{{\tilde{k}}_1}{⟶}\mathrm{P}+\mathrm{E},\hfill \\ \hfill & \mathrm{C}+\mathrm{B}\stackrel{k_2}{⟶}\mathrm{P}+\mathrm{E},\hfill & \hfill & \mathrm{A}+\mathrm{B}+\mathrm{P}\stackrel{{\tilde{k}}_2}{⟶}\mathrm{G}+\mathrm{E}.\hfill \\ \hfill & \mathrm{P}+\mathrm{C}\stackrel{k_3}{⟶}\mathrm{G}.\hfill & \hfill \end{array}$$

The goal is to maximize the profit function $J(\mathit{y},\mathit{u}):=(1143.38x_{P,ss}+25.92x_{E,ss})F_R-76.23F_A-114.34F_B$, where $x_{P,ss},x_{E,ss}$ are the measured steady-state concentrations of products $P,E$, and $F_R$ is the summation of the feed flowrates of reactants A and B, i.e., $F_R:=F_A+F_B$. $F_B$ and the reaction temperature $T_R$ are the manipulated variables ($\mathit{u}=[F_B,T_R]$) with a range $[3,6]\mathrm{k}\mathrm{g} {\mathrm{s}}^{-1}$ and $[343.15,373.15]\mathrm{K}$. The steady-state plant model (${\mathit{F}}_p\left(\mathit{u}\right)$) and the nominal model (${\mathit{F}}_m\left(\mathit{u}\right)$) are:

$${\mathit{F}}_p\left(\mathit{u}\right):=\left(\begin{array}{c}{F}_A-F_R{x}_A-m_R{r}_1\\ F_B-F_R{x}_B-m_R(r_1+r_2)\\ -F_R{x}_C+m_R(2r_1-2r_2-r_3)\\ -F_R{x}_E+2m_R{r}_2\\ -F_R{x}_G+1.5m_R{r}_3\\ -F_R{x}_P+m_R(r_2-0.5r_3)\\ F_R-F_A-F_B\end{array}\right)=0,$$

(11)

and

$${\mathit{F}}_m\left(\mathit{u}\right):=\left(\begin{array}{c}{F}_A-F{x}_A-m_R({\tilde{r}}_1+{\tilde{r}}_2)\\ F_B-F{x}_B-m_R(2{\tilde{r}}_1+{\tilde{r}}_2)\\ -F_R{x}_E+2m_R{\tilde{r}}_1\\ -F_R{x}_G+3m_R{\tilde{r}}_2\\ -F_R{x}_P+m_R({\tilde{r}}_1-{\tilde{r}}_2)\\ F_R-F_A-F_B\end{array}\right)=0,$$

(12)

where $x_i$ is the concentration of component i and $m_R:=$ 2105 kg is the mass of the reaction mixture in the reactor. $\mathit{r}:=\{r_1,r_2,r_3\},\tilde{\mathit{r}}:=\{{\tilde{r}}_1,{\tilde{r}}_2\}$ represent the reaction rates in the plant, and in the nominal model. The reaction rates and their kinetic parameters are defined as:

$$\begin{array}{cccc}\hfill & r_1=k_1{x}_A{x}_B,\hfill & \hfill & k_1=1.660\times {10}^6{\mathrm{e}}^{-6666.7/({\mathrm{T}}_{\mathrm{R}}+273.15)},\hfill \\ \hfill & r_2=k_2{x}_B{x}_C,\hfill & \hfill & k_2=7.212\times {10}^8{\mathrm{e}}^{-8333.3/({\mathrm{T}}_{\mathrm{R}}+273.15)},\hfill \\ \hfill & r_3=k_3{x}_C{x}_P,\hfill & \hfill & k_3=2.675\times {10}^6{\mathrm{e}}^{-\mathrm{11,111}/({\mathrm{T}}_{\mathrm{R}}+273.15)},\hfill \\ \hfill & {\tilde{r}}_1={\tilde{k}}_1{x}_A{x}_B^2,\hfill & \hfill & {\tilde{k}}_1=2.19\times {10}^8{\mathrm{e}}^{-8077.6/({\mathrm{T}}_{\mathrm{R}}+273.15)},\hfill \\ \hfill & {\tilde{r}}_2={\tilde{k}}_2{x}_A{x}_B{x}_P,\hfill & \hfill & {\tilde{k}}_2=4.31\times {10}^{13}{\mathrm{e}}^{-\mathrm{12,438.5}/({\mathrm{T}}_{\mathrm{R}}+273.15)}.\hfill \end{array}$$

For the purpose of illustration of the performance of the proposed proactive perturbation scheme, we consider that ${\tau }_p,{\tau }_d$ and ${\tau }_m$ are $3,12$ and 3 $\min$.

5.2. Simulation Results

The modifier adaptation scheme MAWQA with GMA [26] is used in the simulation study to compare the methods active perturbation around the current input, active perturbation around an estimate of the next input, and the proposed proactive perturbation scheme using the tuning parameters listed in

Table 1. Tuning parameters for simulation study.

No.	Parameter	Williams–Otto	Lithiation Reaction
1.	trust-region scaling $\left(\gamma \right)$	3	3
2.	radius of inner circle $\left(\Delta u\right)$	0.1	0.125
3.	step length for finite differences $\left(\Delta h_u\right)$	0.1	0.2
4.	conditionality limit $\left(\delta \right)$	0.2	0.25
5.	threshold to stop active perturbations $\left(\xi \right)$	0.05	0.1
6.	threshold for unsuccessful iteration $\left(\phi \right)$	0.005	0.005
7.	number of times an unsuccessful input	3	3
	is probed $\left({\chi }_{count}\right)$
8.	number of inputs looked at for	5	5
	convergence $\left(N_c\right)$
9.	accepted variance of J	$0.01$	−
	for convergence $\left(\epsilon \right)$

. In MAWQA with GMA, the plant gradients are estimated from past inputs and steady-state measurements using QA. All schemes are initialized at ${\mathit{u}}^0:=[3,100]$. For the simulation study, the measured plant profit is corrupted by Gaussian noise with zero mean and $1.0$ standard deviation.

Figure 7 illustrates the evolution of the inputs (scaled) and the plant profit function from MAWQA with GMA, MAWQA with GMA with active perturbation schemes APCI and APENI, and MAWQA with GMA with the proposed proactive perturbation scheme. Descriptions of all markers in the figure are given in

Table 2. Description of markers in Figure 7. Input markers in red and black color refer to the input variable $u_1$; blue and magenta input markers refer to the input variable $u_2$.

Marker	Description
[ , ]	successful-iteration input ($[u_1,u_2]$)
[ , ]	unsuccessful-iteration input ($[u_1,u_2]$)
[ , ]	perturbation input ($[u_1,u_2]$)
	value of the plant profit function

. The axis for the value of the profit function ( ) can be found on the right-side of the plot and the axis for the inputs (scaled) to the plant is on the left. Each vertical grid line in the plot represents the completion of an iteration, which includes all steps from giving input to a process to computing a new input by solving a modifier-adaptation problem. For MAWQA with GMA and MAWQA with GMA with active perturbation schemes APCI and APENI, after applying ${\mathit{u}}^0$ ([ , ]) at $t=0$, two input steps (perturbations [ , ]) with a scaled step length $\Delta h_u$ are performed to compute the plant gradients using finite differences. Each input is applied to the plant until it reaches a steady-state, i.e., for ${\tau }_p:=3$ $\min$. The steady-state plant measurement ( ) for each applied input are obtained after ${\tau }_t:=18$ $\min$ from applying an input. Upon computing the plant gradients using the steady-state plant measurements, the modifier adaptation problem in (4) is solved to compute the the input ${\mathit{u}}^1$. As the condition for cardinality is not met, i.e., ${\mathbb{U}}^1<\ell$, $n_u:=2$ additional perturbations ([ , ]) are performed after ${\mathit{u}}^1$ ([ , ]) to compute the plant gradients. Similar to the 0th iteration, the MA problem in (4) is solved in the $1\mathrm{st}$ iteration to compute ${\mathit{u}}^2$. From here on, the cardinality condition is always satisfied. Therefore, QA is used for gradient approximation for all further iterations. In MAWQA with GMA, additional input perturbations are made only if the criticality-check criterion is not satisfied, for example in the $2\mathrm{nd}$ iteration after applying ${\mathit{u}}^2$ to the plant.

In MAWQA with GMA with active perturbation schemes APCI and APENI, additional perturbations are made once the cardinality condition is satisfied, i.e., from the 2nd iteration onwards, and if there are no past inputs lying closer than $\xi$ to the scaled planned perturbation inputs. In MAWQA with GMA with APCI, the additional perturbations are stopped from the 11th iteration, as from $k:=11$, there are no past inputs lying farther than $\xi$ from all the planned perturbation inputs. In MAWQA with GMA with APENI, the perturbation inputs in kth iteration are chosen around an estimate of iteration input for $k+1$th iteration. From the iteration $k:=10$ onwards, all the planned perturbation inputs are not farther than $\xi$ from earlier inputs, so all active perturbations are stopped.

In the new scheme, MAWQA with GMA with proactive perturbation, after applying the input ${\mathit{u}}^0$ ([ , ]) to the plant, as $\varphi :=1$ and $\varphi <\ell$, two input perturbations ([ , ]) ${\mathit{u}}_{f{d}_1}^0,{\mathit{u}}_{f{d}_2}^0$ are performed to compute the plant gradients using finite differences. All necessary steady-state measurements for gradient correction, i.e., ${\mathit{u}}^0$ and two additional perturbations ${\mathit{u}}_{f{d}_1}^0,{\mathit{u}}_{f{d}_2}^0$ are available after $\zeta :=$ 24 $\min$. During the waiting period, i.e., from $(n_u+1)\times max\left\{{\tau }_p,{\tau }_m\right\}:=$ 9 $\min$ to $\zeta :=$ 24 $\min$, $P_{max}:=5$ additional input perturbations ${\mathit{u}}_{a{p}_1}^0,\dots ,{\mathit{u}}_{a{p}_5}^0$ are performed to gain additional plant information. As the input perturbations when $\varphi \ge \ell$ are considered as inputs from MA (according to the flowsheet in Figure 6), among the 5 additional perturbations, ${\mathit{u}}_{a{p}_3}^0,\dots ,{\mathit{u}}_{a{p}_5}^0$ are renamed as ${\mathit{u}}^1,{\mathit{u}}^2,{\mathit{u}}^3$. At 24 $\min$, steady-state measurements for ${\mathit{u}}^0$ and ${\mathit{u}}_{f{d}_1}^0,{\mathit{u}}_{f{d}_2}^0$ are available. Therefore, an MA problem is formulated and solved to compute the next input ${\mathit{u}}^4$. For the next iteration, as $\varphi :=9$ and $\varphi \nless \ell$, FD is no more used for gradient approximation. However, there are not enough measurements for QA, i.e., $\psi :=4$, after applying ${\mathit{u}}^4$. As $(\ell -\psi ):=2$, input perturbations ${\mathit{u}}^5$ and ${\mathit{u}}^6$ are added to smoothly switch from using FD for gradient approximation to QA. After ${\mathit{u}}^6$, $\psi =\ell$ and a new steady-state measurement is always available every $max\left\{{\tau }_p,{\tau }_m\right\}$ thereby overcoming ${\tau }_d$ and ${\tau }_m$. From here on, a new iteration input is computed every $max\left\{{\tau }_p,{\tau }_m\right\}:=3$ $\min$ taking into account the latest measurement information available.

Although all the schemes converged to an input in the neighborhood of the plant optimum, the MAWQA with GMA scheme took the most time, 375 $\min$ and 20 iterations to converge. The time of convergence is indicated with a red-dashed vertical line in Figure 7. The MAWQA with GMA, APCI, and APENI gained additional plant information from suboptimal input perturbations, and therefore converged in 303, 249 $\min$ and $16,13$ iterations. Finally, MAWQA with GMA and proactive perturbation scheme converged in only 123 $\min$ with 31 iterations. It gained more information per time from additional perturbations than that of the earlier proposed active perturbation schemes and also by using the measurement information immediately after it is available.

The mean profit function for 100 simulations for all schemes is shown in Figure 8. The mean and standard deviation of the profit function did not vary significantly upon convergence. The mean and standard deviation of the times of convergence for all modifier adaptation schemes is also shown. In the proactive perturbation scheme, the profit function converges earlier to the plant optimum, followed by the active perturbations schemes APCI, APENI, and the standard modifier adaptation scheme. The mean time of convergence for the MAWQA with GMA if ${\tau }_d=0$, shown in green, is longer than the mean time of convergence of the proposed proactive perturbation scheme as the proposed scheme also overcomes ${\tau }_m$ in addition to ${\tau }_d$. This illustrates that the proposed scheme overcame the effect of time delay $({\tau }_d>>{\tau }_p,{\tau }_m)$ in this case.

The Friedman ranking test [40] was performed for a ranking comparison of the proposed proactive perturbation scheme with standard MAWQA, and the APCI and APENI schemes based on the time of convergence for 100 realizations of the measurement noise. The Friedman test is a nonparametric statistical test to recognize performance differences among different methods across multiple test samples. For each realization of the measurement noise, the methods are ranked based on their time of convergence. The lowest mean rank for a method implies its consistently superior performance, and vice versa for the highest mean rank. The mean rank of the standard scheme with no additional perturbations is $3.87$, for APCI it is $2.6$, for APENI $2.51$, and the proposed proactive perturbation scheme has a mean rank of $1.02$. Thus, the proposed scheme clearly showed a superior performance when compared to that of the other methods.

6. Lithiation Process Case Study

6.1. Process Description

The lithiation case study is a real case where ${\tau }_d>>\{{\tau }_p,{\tau }_m\}$ [41]. The iterative real-time optimization method MAWQA was successfully applied to the real process by Gottu Mukkula et al. [41] but there is still scope for improvement in the time of convergence to the real process optimum due to the large $\mathit{waiting} \mathit{period}$ in the process. In this simulation study, we aim at proving that the plant optimum could be identified much earlier with the proposed proactive perturbation scheme.

A lithiation reaction [41] is performed to produce Lithium 2-Nitrodiphenylamine (component H) in a containerized reactor module developed within the ${\mathrm{F}}^3$-factory project [42] at the INVITE GmbH facility. All chemical components involved in the lithiation reaction are listed in

Table 3. List of chemical components involved in lithiation reaction.

Label	Chemical Name (Abbreviation)	Role
A	Aniline (An)	reactant
B	Lithium bis(trimethylsilyl)amide (Li-HMDS)	reactant
C	Lithium phenylazanide (Li-An)	intermediate
D	Hexamethyldisilazane (HMDS)	by-product
E	1-Fluoro-2-nitrobenzene (FNB)	reactant
F	2-Nitrodiphenylamine (NDPA)	intermediate
G	Lithium fluoride (LiF)	by-product
H	Lithium 2-Nitrodiphenylamine (Li-NDPA)	product
	Tetrahydrofuran (THF)	solvent

and the reaction mechanism is shown in Figure 9. The main modules of the containerized reactor, i.e., the coiled tubular reactor, a NIR flow cell, an online-NMR sensor [43] and a product filter, are shown in Figure 10. The byproduct Lithium fluoride (component G) precipitates along the length of the tubular reactor, and the solid G is filtered out before the reaction mixture from the outlet of the coiled reactor is collected in a product tank. A filtered sample of the reaction mixture is continuously fed to an online-NMR measurement system and to a NIR flow cell to measure the concentration of the product in the reaction mixture at the reactor outlet [43].

The following steady-state plant model was developed assuming a negligible value of the dispersion coefficients for concentrations and temperature:

$$\begin{array}{cc}\hfill \frac{\partial C_i}{\partial z}& =\frac{r_i}{v_z}\forall i\in \{A,B,C,D,E,F,G,H\},\hfill \end{array}$$

(13a)

$$\begin{array}{ll}\frac{\partial T_R}{\partial z}& =\frac{1}{v_z\rho c_p}(KA(T_R-T_{Envi})-{\Delta }_{H_1}{k}_1{C}_A{C}_B-\\ \hfill & {\Delta }_{H_2}{k}_2{C}_C{C}_E-{\Delta }_{H_3}{k}_3{C}_C{C}_F-\\ & {\Delta }_{H_4}{k}_4{C}_B{C}_F),\end{array}$$

(13b)

where $r_i$ is the reaction rate of component i, described by:

$$\begin{array}{cccc}\hfill r_A& =-k_1{C}_A{C}_B+k_3{C}_C{C}_F,\hfill & \hfill & r_E=-k_2{C}_C{C}_E,\hfill \end{array}$$

(14a)

$$\begin{array}{cccc}\hfill r_B& =-k_1{C}_A{C}_B+k_4{C}_B{C}_F,\hfill & \hfill & r_F=k_2{C}_C{C}_E-k_3{C}_C{C}_F-k_4{C}_B{C}_F,\hfill \end{array}$$

(14b)

$$\begin{array}{cccc}\hfill r_C& =k_1{C}_A{C}_B-k_2{C}_C{C}_E-k_3{C}_C{C}_F,\hfill & \hfill & r_G=k_2{C}_C{C}_E,\hfill \end{array}$$

(14c)

$$\begin{array}{cccc}\hfill r_D& =k_1{C}_A{C}_B+k_4{C}_F{C}_B,\hfill & \hfill & r_H=k_3{C}_C{C}_F+k_4{C}_F{C}_B.\hfill \end{array}$$

(14d)

The model parameters are provided in

Table 4. Parameter values used in plant model.

Parameter	Value	Parameter	Value
$\Delta H_1,\dots ,\Delta H_4$	$-6$ $\mathrm{J}$/$\mathrm{mol}$	$T_{in}$	$283.15$ $\mathrm{K}$
$k_{10},\dots ,k_{40}$	$17.8$${\mathrm{s}}^{-1}$	E	$25\times {10}^3$ $\mathrm{J}$/$\mathrm{mol}$
$KA$	$0.001$ $\mathrm{W}$/$\mathrm{K}$	$\Delta H$	$-19$ $\mathrm{J}$/$\mathrm{mol}$
$T_{Envi}$	$293.15$ $\mathrm{K}$	$c_p$	$0.123$ $\mathrm{J}$/($\mathrm{mol}$$\mathrm{K}$)
$\rho$	900 $\mathrm{k}\mathrm{g}$/${\mathrm{m}}^3$	R	$8.314$ $\mathrm{J}$/($\mathrm{mol}$$\mathrm{K}$)

, and the reaction rates constants $k_i\forall i=1,\dots ,4$ depend on the temperature via an Arrhenius-type equation:

$$k_i=k_{i0}exp(\frac{-E}{R{T}_R}).$$

(15)

Although a four-step reaction takes place in the reactor, the nominal model instead considers a single reaction:

$$\mathrm{A}+2\mathrm{B}+\mathrm{E}⟶2\mathrm{D}+\mathrm{G}+\mathrm{H},$$

(16)

leading to a structural plant-model mismatch. The steady-state nominal model is defined as follows:

$$\begin{array}{cc}\hfill \frac{\partial C_i}{\partial z}& =\frac{{\overline{r}}_i}{v_z}\forall i\in \{A,B,D,E,G,H\},\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill \frac{\partial T_R}{\partial z}& =\frac{1}{v_z\rho c_p}\left(KA(T_R-T_{Envi})-{\Delta }_{H}k{C}_A{C}_B{C}_E\right).\hfill \end{array}$$

(17b)

The reaction rates ${\overline{r}}_i$ in the nominal model follow the elementary rate law with the rate constant k defined by the Arrhenius equation with $k=17.8exp(\frac{E}{R{T}_R})$.

The goal is to maximize the plant profit that is computed by the function

$$J(\mathit{y},\mathit{u}):=\frac{{\overline{w}}_H{C}_{H,ss}{M}_H}{\rho }\left(8+\sum _{i=1}^2{u}_i\right)-{\overline{w}}_A{u}_A-{\overline{w}}_E{u}_E,$$

(18)

where $u_A$ and $u_E$ are the feed flowrates of Aniline and FNB in $\mathrm{k}\mathrm{g}$ ${\mathrm{h}}^{-1}$. The feed flow rate of LiHMDS is set to the maximum of 8 $\mathrm{k}\mathrm{g}$ ${\mathrm{h}}^{-1}$. The manipulated variables $u_A$ and $u_E$ can take values between 3 and 8 $\mathrm{k}\mathrm{g}$ ${\mathrm{h}}^{-1}$. $M_H=$ $0.1752$ $\mathrm{k}\mathrm{g}$ ${\mathrm{mol}}^{-1}$ is the molecular weight of the component H, $\{{\overline{w}}_H,{\overline{w}}_A,{\overline{w}}_E\}$ reflect the component prices, their values are {450,000, 10,000, 12,000} $ ${\mathrm{kg}}^{-1}$ and $C_{H,ss}$ is the measured steady state concentration of component H in $\mathrm{mol}$ ${\mathrm{m}}^{-3}$. $\{{\tau }_p,{\tau }_d,{\tau }_m\}$ are $\{3,12,3\}$ $\min$.

6.2. Simulation Results

All tuning parameters used in the simulation study are listed in Table 1. The RTO schemes are initialized at ${\mathit{u}}^0:=[3.58,3.58]$ and the measured variable $C_H$ is corrupted by Gaussian noise with zero mean and $1.0$ standard deviation. Similar to the Williams–Otto reactor case study, model adequacy is enforced using the guaranteed model adequacy (GMA) proposed in Gottu Mukkula and Engell [26]. Convergence to the process optimum, in this case study, is assumed when the profit from $N_c$ consecutive iterations is greater than $1.175\times {10}^5 \$$. Figure 11 illustrates the mean and the 10 to 90 percentile range of the profit function, and the mean time of convergence for 100 simulations of MAWQA with GMA, MAWQA with GMA with active perturbation schemes APCI and APENI, and MAWQA with GMA with the proposed proactive perturbation scheme. The deviations in the profit function are due to the measurement noise, which affects the computation of the next inputs via the calculation of the quadratic approximation. The figure shows that the proposed proactive perturbation scheme converges faster when compared to that of MAWQA with GMA and the two active perturbation variants proposed earlier. The reason for this is that the idle time is used better by intelligently scheduling the process inputs taking into account the time delays to obtain the steady-state process measurements. Figure 12 shows the mean value of the profit function from MAWQA with GMA, MAWQA with GMA, APCI, and APENI, and MAWQA with GMA and the proposed proactive perturbation scheme. The mean time of convergence and the standard deviation of the time of convergence for 100 simulations are shown as error bars. In addition to the faster convergence of the proposed proactive perturbation scheme, the standard deviation is also smaller when compared to the other schemes. Similar to the Williams–Otto reactor case study, the mean time of convergence of the proposed proactive perturbation scheme is even shorter than that of the MAWQA with GMA scheme with ${\tau }_d=0$. The mean Friedman rank [40] of the standard scheme with no additional perturbations is $3.64$, for APCI it is $2.82$, for APENI $2.5$, and for the proposed proactive perturbation scheme it is $1.06$.

The time to compute a new input in an Intel(R) Core(TM) i7-4790 CPU @ 3.60 GHz × 4 using a 64 bit Windows 10 operating system for both case studies ranged between 0.01 and 0.02 s. Such short computation time would enable the usage of the proposed proactive perturbation scheme also for processes with very short settling time.

7. Conclusions

In this paper, we address iterative RTO with plant-model mismatch in the presence of a time delay caused due to the positioning of an analytic measurement device at some distance from the place where the process stream is sampled. The time delay leads to a longer time of convergence to the real process optimum for modifier adaptation methods, during which the plant is not operated optimally. Perturbation schemes to handle this time delay in the sense that the period where one must wait for the result of the transfer of the sample and the analytic measurement is used for further probing of the plant were proposed earlier, but they rely on heuristics for the choice of the plant perturbations. A new proactive perturbation scheme is proposed in this paper to handle the time delay where additional perturbation inputs are computed by solving an optimization problem using all available process information, which contrasts the earlier proposed active perturbation schemes.

In the proposed proactive perturbation approach, the measurement delay affects the time of convergence only once at the beginning, until a number of plant measurements for gradient approximation using quadratic approximation are available. Later, the proactive perturbation scheme computes a new input by solving a modifier adaptation problem as soon as new measurement information is available, with a sampling time equal to the maximum of the settling time of the plant and the time taken by the measurement device to analyze a sample. The performance of the proposed scheme was evaluated for two examples and compared with that of other active perturbation schemes using two case studies. The proposed proactive perturbation scheme drives the real plant to its optimum significantly faster than the other approaches, and it significantly reduces the effect of the measurement time delay for both case studies because the scheme collects more useful information from the perturbations. An important issue in the application of MA schemes to real plants is to freeze and restart the adaptation when the input moves are mainly caused by measurement noise. Initial proposals for this can be found in [38], but this requires more research in the future.