Dynamic and Statistical Operability of an Experimental Batch Process

The operability approach has been traditionally applied to measure the ability of a continuous process to achieve desired specifications, given physical or design restrictions and considering expected disturbances at steady state. This paper introduces a novel dynamic operability analysis for batch processes based on classical operability concepts. In this analysis, all sets and statistical region delimitations are quantified using mathematical operations involving polytopes at every time step. A statistical operability analysis centered on multivariate correlations is employed for the first time to evaluate desired output sets during transition that serve as references to be followed to achieve the final process specifications. A dynamic design space for a batch process is, thus, generated through this analysis process and can be used in practice to guide process operation. A probabilistic expected disturbance set is also introduced, whereby the disturbances are described by pseudorandom variables and disturbance scenarios other than worst-case scenarios are considered, as is done in traditional operability methods. A case study corresponding to a pilot batch unit is used to illustrate the developed methods and to build a process digital twin to generate large datasets by running an automated digital experimentation strategy. As the primary data source of the analysis is built in a time-series database, the developed framework can be fully integrated into a plant information management system (PIMS) and an Industry 4.0 infrastructure.


Introduction
The concepts of process operability defined by Vinson and Georgakis [1] have been widely studied and applied to several linear and non-linear systems. In essence, process operability allows the analysis of the achievable output set (AOS) that a process can reach, considering an available input set (AIS) and an expected disturbance set (EDS), for a given process model. This analysis can also verify if the desired output set (DOS) can be achieved and quantify how much of the AOS covers the DOS specifications by defining an operability index (OI) [2]. The inputs for the analysis correspond to the process design and operational limitations, environmental constraints, as well as a group of external disturbances. The outputs cover production targets, quality specifications, business benchmarks, and safetycritical conditions [3].
Past operability studies have also considered the application of operability to control strategies [4,5], process intensification [6,7], and design optimization [8,9]. The operability framework has been extensively studied from a steady-state point of view [10], in which the graphical regions used consisted of mainly stable points of operation [11,12]. Other studies have investigated the dynamic operability concerning the minimum time to achieve different steady states [13,14]. Additionally, these previously reported studies have considered only deterministic sets, without statistical or probabilistic concepts being involved. The process model (M) can be defined by the first principles or process simulators (e.g., in Aspen Plus, DYNSIM, AVEVA Process Simulation, etc.) to reproduce a real unit and obtain the relationships between input and output spaces [3]. The following sets are essential for the operability analysis: a. Available input set (AIS)-Considers all available inputs that may change within a specific range according to physical or operating constraints. This set may represent design or operational inputs. Operational inputs are manipulated variables (MV) subject to control studies, whereas design inputs can be related to process materials, design dimensions, etc. [11,16]. This set is defined as: b. Expected disturbance set (EDS)-Represents the known process uncertainties, usually chosen as the extreme values in order to consider worst-case scenarios [16].
c. Achievable output set (AOS)-Compiles all possible outputs that the system can achieve, given the AIS and considering an EDS. Thus, this set is defined as a function of inputs and disturbances. For each disturbance scenario, ∈ , the AOS can be obtained as [12]: The intersection of all AOS(d)s yields the AOS (shown in Equation (4)), which is a region where the setpoints of the output variables can be moved successfully using the available control action in the AIS while compensating for any disturbances in the EDS [12]. = ∈ (4) d. Desired output set (DOS)-Corresponds to the desired region to be obtained for the process. This set may be defined, for example, by production requirements such as production rates or product qualities, process specifications, or efficiency. The DOS is mathematically defined in Equation (5) [11]. The process model (M) can be defined by the first principles or process simulators (e.g., in Aspen Plus, DYNSIM, AVEVA Process Simulation, etc.) to reproduce a real unit and obtain the relationships between input and output spaces [3]. The following sets are essential for the operability analysis: a.
Available input set (AIS)-Considers all available inputs that may change within a specific range according to physical or operating constraints. This set may represent design or operational inputs. Operational inputs are manipulated variables (MV) subject to control studies, whereas design inputs can be related to process materials, design dimensions, etc. [11,16]. This set is defined as: Expected disturbance set (EDS)-Represents the known process uncertainties, usually chosen as the extreme values in order to consider worst-case scenarios [16].
c. Achievable output set (AOS)-Compiles all possible outputs that the system can achieve, given the AIS and considering an EDS. Thus, this set is defined as a function of inputs and disturbances. For each disturbance scenario, d ∈ EDS, the AOS can be obtained as [12]: The intersection of all AOS(d)s yields the AOS (shown in Equation (4)), which is a region where the setpoints of the output variables can be moved successfully using the available control action in the AIS while compensating for any disturbances in the EDS [12]. d.
Desired output set (DOS)-Corresponds to the desired region to be obtained for the process. This set may be defined, for example, by production requirements such as production rates or product qualities, process specifications, or efficiency. The DOS is mathematically defined in Equation (5) [11]. The DOS can be described by straight edges that connect the vertices given by y min j and y max j for each output or discontinuous or disjointed points that discretize such edges [3]. In practice, a process is fully operable if all outputs satisfy the desired requirements, i.e., the DOS is a subset of the AOS. The operability index (OI) is defined to quantify the achievability of the DOS, as shown in Equation (6). If some regions of the DOS are not achievable, the OI is, thus, less than one [12].
where µ represents a measure of the size of the set, e.g., area for 2-D sets, hypervolume for high-D sets, etc. Thus, mathematical operations involving intersections of polytopes, such as AOS ∩ DOS, need to be performed to determine the OI. For such operations, the Multi-Parametric Toolbox (MPT) [22] in MATLAB can be used for convex regions. For non-convex regions, the MATLAB built-in subroutine "boundary" may be employed.
In the classical operability methods, these sets are mainly studied at steady-state conditions. The existing dynamic approach investigates the operability from the standpoint of the minimum time required for a system to move between steady states and reject expected disturbances. Existing dynamic operability methods can also be used to calculate the amount of constraint relaxation necessary to prevent the occurrence of transient infeasibilities when disturbances affect the process [13,16]. The analysis proposed in the next section considers the extension of the existing operability methods for the calculation of dynamic sets employing concepts of statistics.

Proposed Approach: Dynamic and Statistical Operability
The classical steady-state operability concepts have inspired the proposed method introduced in this section, which now considers the entire time frame, including the process transient operations. In particular, the AIS remains the same as the concept described by Equation (1), with the values indexed at the initial time (t 0 ). The other sets have to be redefined, considering their categorization in terms of their variable abilities to be manipulated during process operation. For example, a reactant concentration can be manipulated at the beginning of the batch, however once the reaction occurs, it will be dependent on the process design and mass and energy balances. However, the amount of fluid entering the jacket is always changeable. Therefore, two categories may be defined for the process variables: a.
Static-no handling is available during the process and it is typically defined at the beginning, e.g., for the initial temperature and concentrations; b.
Dynamic-the values can be changed during operation, whether in an open or closed loop, e.g., for the thermal fluid inlet flow and temperature.

Expected Disturbance Set
The new EDS definition considers a probabilistic approach, in which not only are the worst-case scenarios analyzed, but also the intermediate ones with more likelihood of occurring. The dynamic and probabilistic expected disturbance set (EDS p (t)) is defined by Equation (7) as follows: • Disturbance intensity-The deviation from the nominal value of the disturbance variable is assumed random for unbiased processes and with a Gaussian behavior. Essentially, this work considers ±3σ limits for d min k and d max k as in Equation (7) and the generation of random disturbances within these bounds;

•
Dynamic expectation-A disturbance may occur at any moment during operations, reaching a maximum value and then reducing its magnitude due to a control action; Processes 2021, 9, 441 5 of 18 thus, it can have different behaviors (e.g., step or impulse). In this work, it is assumed that the disturbance follows a ramp, reaching the maximum value at 25%, 50%, or 75% for a one-hour long batch. Such disturbance intensity then remains for 5, 15, or 30 min, respectively. Figure 2 illustrates a disturbance example in which the temperature is affected, reaching 61.52 • C at the 15 min mark, remaining for 5 min, and then returning to its nominal value of 60.00 • C.
Processes 2021, 9, x FOR PEER REVIEW 5 of 20 • Dynamic expectation-A disturbance may occur at any moment during operations, reaching a maximum value and then reducing its magnitude due to a control action; thus, it can have different behaviors (e.g., step or impulse). In this work, it is assumed that the disturbance follows a ramp, reaching the maximum value at 25%, 50%, or 75% for a one-hour long batch. Such disturbance intensity then remains for 5, 15, or 30 min, respectively. Figure 2 illustrates a disturbance example in which the temperature is affected, reaching 61.52 °C at the 15 min mark, remaining for 5 min, and then returning to its nominal value of 60.00 °C. The more disturbance possibilities that are considered, the more representative the EDSp(t) is. In real batch reaction systems, historical data can be investigated to further characterize the disturbance behavior throughout the operation. In a PIMS infrastructure, time-series data for different variables, including disturbances, may come from sensors, estimators, or historical databases.

Achievable Output Set
As the outputs are considered time-dependent, AOS(d,t) defined by Equation (8) is indexed per timestamp and by the disturbance d, which is part of the EDSp(t).

, =
∈ ℝ | = , , ; for all ∈ and a fixed ∈ The AOS(t) is similarly calculated as in the original definition, based on the intersection of the AOS(d,t)s for a specific time t given by Equation (9).

Desired Output Set
Previous operability studies shows the desired output set (DOS) as a set desired to be achieved at steady state [3,9]. For intrinsic time-dependent batch processes, the interest The more disturbance possibilities that are considered, the more representative the EDS p (t) is. In real batch reaction systems, historical data can be investigated to further characterize the disturbance behavior throughout the operation. In a PIMS infrastructure, time-series data for different variables, including disturbances, may come from sensors, estimators, or historical databases.

Achievable Output Set
As the outputs are considered time-dependent, AOS(d,t) defined by Equation (8) is indexed per timestamp and by the disturbance d, which is part of the EDS p (t).
for all u ∈ AIS and a fixed d ∈ EDS p (t) (8) The AOS(t) is similarly calculated as in the original definition, based on the intersection of the AOS(d,t)s for a specific time t given by Equation (9).

Desired Output Set
Previous operability studies shows the desired output set (DOS) as a set desired to be achieved at steady state [3,9]. For intrinsic time-dependent batch processes, the interest lies in the final process specification. Thus, the DOS(t f ) can be defined, where t f denotes the final time for the same definition as in Equation (5).
The proposed dynamic operability framework has to first consider the conditions capable of satisfying the DOS(t f ). Based on these conditions, the dynamic profiles can be mapped backwards and the desired output set can be evaluated per timestamp, defining the DOS(t). This mapping strategy is based on a joint reliability region, which carries information about all variables/dimensions of DOS(t f ) and their self-correlations based on the covariance matrix. The proposed method, thus, relies on multivariate distributions that change their shape according to the variable interactions. Figure 3 is based on the concepts presented in [18,23] and shows an example from a bidimensional analysis of the desired regions represented by an ellipse.
Processes 2021, 9, x FOR PEER REVIEW 6 of 20 lies in the final process specification. Thus, the DOS( ) can be defined, where denotes the final time for the same definition as in Equation (5).
The proposed dynamic operability framework has to first consider the conditions capable of satisfying the DOS( ). Based on these conditions, the dynamic profiles can be mapped backwards and the desired output set can be evaluated per timestamp, defining the DOS(t). This mapping strategy is based on a joint reliability region, which carries information about all variables/dimensions of DOS( ) and their self-correlations based on the covariance matrix. The proposed method, thus, relies on multivariate distributions that change their shape according to the variable interactions. Figure 3 is based on the concepts presented in [18,23] and shows an example from a bidimensional analysis of the desired regions represented by an ellipse. For the ellipse in Figure 3, a perfect circle will be found if there is no correlation between and , while the ellipse becomes "thinner" as the relationship between variables strengthens. This shape change in Figure 3 is due to the eigenvalues ( ) present in the semi-axis expression, which can be evaluated by Equation (10).
in which , is the critical value of a chi-square distribution with degrees of freedom and % level of significance. Additionally, the shape geometry is centered around the mean value ( , ) and has directions defined by the eigenvectors ( and ) and the angle . The ellipse parametric coordinates are described by Equations (11) and (12).
These multivariate concepts are the basis of the DOSp(t) definition. The DOSp(t), thus, represents the desired output region at a specific time t, in which the "p" subscript is due to its probabilistic nature. Therefore, the DOSp(t) can be mathematically defined as in Equation (13). For the ellipse in Figure 3, a perfect circle will be found if there is no correlation between y 1 and y 2 , while the ellipse becomes "thinner" as the relationship between variables strengthens. This shape change in Figure 3 is due to the eigenvalues (λ j ) present in the semi-axis expression, which can be evaluated by Equation (10).
in which χ 2 ν,α is the critical value of a chi-square distribution with ν degrees of freedom and α% level of significance. Additionally, the shape geometry is centered around the mean value (y 1 , y 2 ) and has directions defined by the eigenvectors (e 1 and e 2 ) and the angle θ. The ellipse parametric coordinates are described by Equations (11) and (12).
These multivariate concepts are the basis of the DOS p (t) definition. The DOS p (t), thus, represents the desired output region at a specific time t, in which the "p" subscript is due to its probabilistic nature. Therefore, the DOS p (t) can be mathematically defined as in Equation (13). where n is the DOS dimension; thus, when n = 2, Equation (13) corresponds to an ellipse, n = 3 is an ellipsoid, and so on. Finally, the dynamic operability index (OI(t)) is evaluated at each timestamp according to Equation (14). The lower the OI(t), the more attention is needed in the specific analyzed time period to guarantee the desired process specification. If no process condition can satisfy the DOS(t f ), the OI(t f ) is null, resulting in batch run conditions that are not able to achieve the desired specifications.

Case Study Definition and Setup
This section presents a case study that consists of a pilot batch reactor system. The proposed operability script for this system is described along with the time-series dataset used. Additionally, the developed dynamic model for this system is validated against the pilot bench unit, producing a digital twin for the process. Following the script, automatic model runs are performed considering all of the recipes within a digital experimentation strategy encompassing 900 simulated conditions (with approximately 1000 generated points per simulated condition), which are defined based on the AIS and the EDS p (t).

Pilot Batch Unit
The addressed system is characterized by an endothermic second-order reaction between sodium hydroxide (NaOH) and ethyl acetate (CH 3 COOCH 2 CH 3 ), as described in Equation (15). This reaction is carried out in a 1.2L stirred batch reactor and follows the kinetic parameters according to [24]. In this reactor system, water flows in a jacket to transfer enough heat for isothermal operation. A pH indicator is used to mark the end of the reaction process. The experimental apparatus is depicted in Figure 4, where a thermocouple and an electrical conductivity meter are used to acquire the reactor's temperature and concentration, respectively.
Processes 2021, 9, x FOR PEER REVIEW 8 of 20 The heating fluid flows in a closed circuit, where electrical resistance is manipulated to control the reactor temperature. A peristaltic pump guarantees its setpoint flow, which is seated in a supervisory control and data acquisition (SCADA) system present in the vendor software and used in the pilot unit operation. The system is also able to acquire  The heating fluid flows in a closed circuit, where electrical resistance is manipulated to control the reactor temperature. A peristaltic pump guarantees its setpoint flow, which is seated in a supervisory control and data acquisition (SCADA) system present in the vendor software and used in the pilot unit operation. The system is also able to acquire the data as a time-series database per batch. The reactants have been previously prepared at determined concentrations and volumes per experiment. The conversion is calculated based on the conductivity measurements, and along with temperature, represents the most relevant variables in the study.
The digital experimentation strategy is delimited by a sensitivity analysis involving different ethyl acetate concentrations (0.1, 0.05, and 0.07 mol/dm 3 ) and temperatures (30, 35, and 40 • C), with a fixed sodium hydroxide concentration of 0.1 mol/dm 3 . The experimental pilot plant was used to develop and validate a process digital twin, enabling automated and extensive digital manipulation of the system with minimal financial cost and safety risks.

Dynamic Operability Sets Applied to the Batch Case Study
The categorization of the process variables employed is detailed in Table 1, and was used to define the available input set (AIS), as well as the probabilistic expected disturbance set (EDS p (t)). This categorization is essential when studying the AIS and EDS.

Category Variable Description
Static u 1 Ethyl acetate initial molar concentration (C EtAc,0 ) u 2 Sodium hydroxide initial molar concentration (C NaOH,0 ) u 3 Initial reactor temperature (T 0 ) d 1 Sodium hydroxide initial volume (V NaOH,0 ) The considered AIS has four dimensions, as given in Equation (16), with three static variables, i.e., ones that can only be defined at the beginning of the batch. In particular, the two reactant solutions were previously prepared at specified concentrations. Additionally, the initial reactor temperature depends on the environmental conditions and the heating system capacity.
Equation (17) describes the EDS p . Unlike d 3 , the two static disturbances (d 1 and d 2 ) only need an intensity to be defined in the beginning of the batch, with no duration. A ramp behavior is assumed for the dynamic disturbance to add complexity and simulate the real case, in which the disturbance first changes from its nominal value and varies to a maximum, where it remains for a certain period until a control system brings it back to its nominal value. The specifically considered disturbance values are detailed in Table A1 in the Appendix A.
Processes 2021, 9, 441 9 of 18 The DOS(t f ) is defined by Equation (18), where y 1 and y 2 correspond to reactor conversion of sodium hydroxide and temperature, respectively. This definition is similar to Equation (5), but related to the batch end time. The batches that satisfy the DOS(t f )-also known as the deterministic desired output set in this work-will serve as basis to calculate the intermediate DOS p (t)s, i.e., the probabilistic DOS's.

Automated Digital Experimentation
Due to the large number of experiments to be performed, an automated digital experimentation strategy is proposed using the Aspen software capability to communicate with Excel and Visual Basic for Applications (VBA) via Aspen Simulation Workbook (ASW). The batch recipes are defined by a multivariable experimental design formed by combinations of AIS and EDS p (t), with variables and conditions tabulated in a database. The strategy logics and queue organization are described in the experimentation planning (EP) module in Figure 5, in which a VBA code is written to assemble the recipe specifications and pass them to the process digital twin (DT) via ASW ( Figure 5-step I). The MATLAB code (described in the next section) and VBA dynamic arrays organize the data in a tridimensional way according to Figure 6, where k response variables along time to have the information for the n batches. In such a way, all developed scripts can be expanded to more extensive cases as necessary. The DT corresponds to an Aspen (ABM) model that provides the results at different timestamps because the AspenTech internal numerical method adjusts the time steps per batch to provide better convergence; hence, the data need to be reconciled before being archived in the final dataset. For this study, around 100 points for 10 variables are stored per one-hour batch. The outputs are sent to the database manager (DM) via ASW ( Figure 5-step II). The DM is another VBA script responsible for writing the simulation outputs that converge in a time-series database. Not necessarily all runs will converge, since the EDS p (t) is pseudo-randomly generated; thus, a condition with no physical solution could be found, which would be reported with a trial status. Once the batch reaches the end, a signal is sent to EP via ASW ( Figure 5-step III), and then the ABM is reset to start the next simulation. This loop detailed in Figure 5 is repeated until the last experiment is performed. The complete automated digital experimentation process produces approximately 900,000 data points.
The MATLAB code (described in the next section) and VBA dynamic arrays organize the data in a tridimensional way according to Figure 6, where k response variables along time t 1 to t f have the information for the n batches. In such a way, all developed scripts can be expanded to more extensive cases as necessary. The MATLAB code (described in the next section) and VBA dynamic arrays organize the data in a tridimensional way according to Figure 6, where k response variables along time to have the information for the n batches. In such a way, all developed scripts can be expanded to more extensive cases as necessary.

Dynamic Operability Algorithm
The dynamic operability analysis is implemented in MATLAB, according to the diagram described in Figure 7. Essentially, the time-series database is loaded first, the variables are identified, and after a few inputs to the algorithm have been defined, the methodology evaluation is started. These inputs correspond to:

Dynamic Operability Algorithm
The dynamic operability analysis is implemented in MATLAB, according to the diagram described in Figure 7. Essentially, the time-series database is loaded first, the variables are identified, and after a few inputs to the algorithm have been defined, the methodology evaluation is started. These inputs correspond to:

Results
Initially, the validation of the digital twin model with the experimental results from the pilot plant is presented. Then, the dynamic operability results are shown for different sets, namely AOS(d,t), DOS(t f ), and DOS p (t), as well as the OI(t) calculations.

Batch Model Validation with Experimental Results
The batch reactor temperature and sodium hydroxide molar concentration are used as variables for the model validation. The experiments are conducted until the acid-base indicator (crystal violet) marks the end of the reaction, with the same experiment duration applied as for the simulation.
The comparisons between the pilot unit data and the simulated digital twin outcomes are presented for the experiments in Figures 8 and 9. The sodium hydroxide concentration and reactor temperature profiles are shown in Figure 8a,b respectively for the fixed initial concentration of ethyl acetate specified in the legend.
The digital twin result profiles in Figure 8 thus closely represent the experimental ones for the entire batch. Specifically, the higher the initial concentration of ethyl acetate (EtAc), the more sodium hydroxide is converted; hence, the sodium hydroxide molar concentration is lowest when ethyl acetate starts at 0.10 mol/L. The highest relative error (4.0% approximately) occurs at the temperature profile for the lowest ethyl acetate initial concentration (0.05 mol/dm 3 ); however, this condition did not substantially affect the molar concentration of sodium hydroxide over time.
as variables for the model validation. The experiments are conducted until the acid-base indicator (crystal violet) marks the end of the reaction, with the same experiment duration applied as for the simulation.
The comparisons between the pilot unit data and the simulated digital twin outcomes are presented for the experiments in Figures 8 and 9. The sodium hydroxide concentration and reactor temperature profiles are shown in Figures 8a and 8b respectively for the fixed initial concentration of ethyl acetate specified in the legend. as variables for the model validation. The experiments are conducted until the acid-base indicator (crystal violet) marks the end of the reaction, with the same experiment duration applied as for the simulation. The comparisons between the pilot unit data and the simulated digital twin outcomes are presented for the experiments in Figures 8 and 9. The sodium hydroxide concentration and reactor temperature profiles are shown in Figures 8a and 8b respectively for the fixed initial concentration of ethyl acetate specified in the legend.  Figures 8 and 9 indicate that the model matches the behavior shown by the experimental results of the pilot unit for all the sensitivity analyses performed for both concentration and temperature profiles. The highest relative error obtained was around 15% when the molar concentration was below 0.02 and the temperature was at 30 • C (see Figure 9), which is likely related to sensor precision. Thus, as the digital twin (DT) parameters are well adjusted to represent the entire batch system, this DT is used in the proposed dynamic and statistical operability approach. This ABM model is then simulated multiple times considering the cases determined by the design experimentation (see summary in Table A1). In this work, nearly a million data points were generated by the DT, with the corresponding data stored in a process database.

Dynamic and Statistical Operability Results
As at t = 0 h, there is no change in the system; the first AOS(d,t) is analyzed at t = 0.1 h, as depicted in Figure 10. In this figure, all of the outputs from the simulations at this timestamp are represented by the markers, while the non-disturbance case, i.e., AOS (d = 0, t = 0.1 h), is represented by the solid black line. The results from the non-disturbance case are plotted along with the ones obtained by the #1 disturbance case from Table A1-generating the AOS (d = 9.7, t = 0.1 h)-in Figure 10a.
(see summary in Table A1). In this work, nearly a million data points were generated by the DT, with the corresponding data stored in a process database.

Dynamic and Statistical Operability Results
As at t = 0 h, there is no change in the system; the first AOS(d,t) is analyzed at t = 0.1 h, as depicted in Figure 10. In this figure, all of the outputs from the simulations at this timestamp are represented by the markers, while the non-disturbance case, i.e., AOS (d = 0, t = 0.1 h), is represented by the solid black line. The results from the non-disturbance case are plotted along with the ones obtained by the #1 disturbance case from Table A1generating the AOS (d = 9.7, t = 0.1 h)-in Figure 10a. Then, the outputs, AOS(d, t = 0.1 h), considering all the disturbance instances (in Table A1) are outlined by the dashed lines (light blue) in Figure 10b. Finally, the AOS(t = 0.1 h) is delimited by the solid line (dark blue) in Figure 10b, where the AOS(t = 0.1 h) is the intersection of all AOS(d,t = 0.1 h)s at 0.1 h and is the one used for operability index (OI) calculation purposes. For this figure, each region was outlined using the MATLAB boundary function, while mathematical operations such as intersections were determined using another built-in MATLAB function called "polybool".
The procedure above is used to determine the AOS(d,t)s presented in Figure 11a with the output points and AOS(t)s in Figure 11b, at all ten timestamps considered during the Then, the outputs, AOS(d, t = 0.1 h), considering all the disturbance instances (in Table A1) are outlined by the dashed lines (light blue) in Figure 10b. Finally, the AOS(t = 0.1 h) is delimited by the solid line (dark blue) in Figure 10b, where the AOS(t = 0.1 h) is the intersection of all AOS(d,t = 0.1 h)s at 0.1 h and is the one used for operability index (OI) calculation purposes. For this figure, each region was outlined using the MATLAB boundary function, while mathematical operations such as intersections were determined using another built-in MATLAB function called "polybool".
The procedure above is used to determine the AOS(d,t)s presented in Figure 11a with the output points and AOS(t)s in Figure 11b, at all ten timestamps considered during the batch operation. The developed script can consider as many intervals as are available in the data; however, for visualization purposes, the analysis here is delimited to 10 time instants. The plot of the deterministic desired output set (DOS( )) at the end of the batch analyzes whether the designed process can achieve the specifications even with the expected disturbances. By deterministic, it is referred here to the set described by Equation (18) The plot of the deterministic desired output set (DOS(t f )) at the end of the batch analyzes whether the designed process can achieve the specifications even with the expected disturbances. By deterministic, it is referred here to the set described by Equation (18), outlined by lower and upper limits and considering no statistical interaction among variables. The output points that satisfy the deterministic DOS(t f ) are indexed and used to build the primary joint reliability region at the final timestamp, i.e., DOS p (t f ), as shown in Figure 12. The plot of the deterministic desired output set (DOS( )) at the end of the batch analyzes whether the designed process can achieve the specifications even with the expected disturbances. By deterministic, it is referred here to the set described by Equation (18), outlined by lower and upper limits and considering no statistical interaction among variables. The output points that satisfy the deterministic DOS( ) are indexed and used to build the primary joint reliability region at the final timestamp, i.e., DOSp( ), as shown in Figure 12. The DOS p (t f ) is, thus, the tolerance region considered for control purposes, with its center pulled by the region with the highest concentration of data, as shown by the ellipse in Figure 12. Additionally, the DOS p (t f ) is outlined by considering the tolerance at 60%. The data in Figure 12 are highly concentrated (with many points overlapping each other) at 98-100% conversion and 58-65 • C, as well as at 35-45 • C and 95-100% conversion, which is the set used to build the DOS p (t f ). Thus, if the desired output set is centered in a range different from these conditions (such as around 47-56 • C and 95-100%), the desired specifications might not be achieved after a 1 h batch with the AIS and EDS p (t) described in this work.
With all outputs indexed by timestamp, the other DOS p (t)s can be calculated from this final point, with the indexed data shown in Figure 13.
In the magnified portion in Figure 13, the shape of the joint reliability region is still closer to an ellipse at the end of the batch, which indicates a strong correlation between temperature and conversion. The funnel-like behavior shown in the figure indicates more desirable possibilities with a higher probability of occurrence at the beginning of the process, followed by the reduction of such possibilities at each time step. This response is due to the fact that the batch is a time-dependent process that relies on cumulative statechanging behavior along the operation (t). Such a process, thus, becomes more difficult to operate as it evolves, as also reflected in the decreasing OI(t) values in Table 2.
The data in Figure 12 are highly concentrated (with many points overlapping each other) at 98-100% conversion and 58-65 °C, as well as at 35-45 °C and 95-100% conversion, which is the set used to build the DOSp( ). Thus, if the desired output set is centered in a range different from these conditions (such as around 47-56 °C and 95-100%), the desired specifications might not be achieved after a 1 h batch with the AIS and EDSp(t) described in this work.
With all outputs indexed by timestamp, the other DOSp(t)s can be calculated from this final point, with the indexed data shown in Figure 13. In the magnified portion in Figure 13, the shape of the joint reliability region is still closer to an ellipse at the end of the batch, which indicates a strong correlation between temperature and conversion. The funnel-like behavior shown in the figure indicates more desirable possibilities with a higher probability of occurrence at the beginning of the process, followed by the reduction of such possibilities at each time step. This response is due to the fact that the batch is a time-dependent process that relies on cumulative statechanging behavior along the operation (t). Such a process, thus, becomes more difficult to operate as it evolves, as also reflected in the decreasing OI(t) values in Table 2.  Similarly to the classical operability analysis, the intersections between the AOS(t)s and DOS p (t)s can be evaluated at each timestamp over the batch campaign to provide a graphical idea (as shown in Figure 14) of how much the desired output set is achievable and where to operate the process in order to guarantee the final specifications would be reachable. In particular, the black marker in Figure 14 exemplifies an operation sample of a batch in which the specifications were achieved, even considering the random disturbances. Thus, this process can be characterized as fully dynamically operable in terms of reactor conversion and temperature specifications. In case of deviation from the dashed area in Figure 14, the series of regions would serve as a path for the operation to drive the process back to the proper DOSp(t). In order to drive the process back to the desired region in this case, a control system (such as classical feedback or advanced model predictive control) would be required.

Conclusions
The dynamic operability approach introduced here was applied to a strictly timedependent process. The developed method consists of an extension of the classical operability concepts to a scenario where process dynamics and disturbance statistics are considered. The statistical concepts explored included multivariate normal distributions and probabilistic analysis of the expected disturbance set (EDS). A new desired output set (DOSp(t)) was proposed based on the joint reliability region in which the tolerance value α definition is necessary. A procedure to generate a dynamic design space for a batch process in the form of a sequence of regions was provided and can be used in practice to assist batch process operations, monitoring, and control. Such regions consist of desirable and achievable sets at each time step, considering the constraints and the disturbances being assumed as normally distributed. The product specifications are guaranteed at the end of the process if operated inside the intersection series between the AOS and DOS regions. Figure 14. Dynamic operability analysis: AOS(t)s, DOS p (t)s, and their intersections over the batch operation. The shaded areas along the campaign represent how much the desired set is achievable.
In case of deviation from the dashed area in Figure 14, the series of regions would serve as a path for the operation to drive the process back to the proper DOS p (t). In order to drive the process back to the desired region in this case, a control system (such as classical feedback or advanced model predictive control) would be required.

Conclusions
The dynamic operability approach introduced here was applied to a strictly timedependent process. The developed method consists of an extension of the classical operability concepts to a scenario where process dynamics and disturbance statistics are considered. The statistical concepts explored included multivariate normal distributions and probabilistic analysis of the expected disturbance set (EDS). A new desired output set (DOS p (t)) was proposed based on the joint reliability region in which the tolerance value α definition is necessary. A procedure to generate a dynamic design space for a batch process in the form of a sequence of regions was provided and can be used in practice to assist batch process operations, monitoring, and control. Such regions consist of desirable and achievable sets at each time step, considering the constraints and the disturbances being assumed as normally distributed. The product specifications are guaranteed at the end of the process if operated inside the intersection series between the AOS and DOS regions.
The process model developed for the performed batch case study was validated with experiments from a pilot plant. This model was defined as the digital twin for the batch process, and an automated digital experimentation strategy was proposed in an Industry 4.0 framework. The success of the analysis performed here depends on the data source reliability and quality generated by the digital twin. The methodology accuracy is improved as more data are available, and from this perspective a plant information management system (PIMS) can be a powerful tool for operability studies. In future work, the authors plan to integrate the current solution with alarms and message notifications whenever the dynamic operability methodology detects operational actions necessary for the process to return to an operable state.