**2. Methodology**

*2.1. Description of the Mathematical Model*

The dynamic model of the segmented fluidised bed dryer being explored here is implemented in the gPROMS modelling suite as part of the gPROMS FormulatedProducts® library [31]. The underlying mathematical formulation is based on the mechanistic model presented by Burgschweiger et al. [35,36] and model parameters have been validated using the Diamond Pilot Plant (DiPP) at the University of Sheffield. For the sake of brevity, we omit the presentation of the full mathematical model and the interested reader is referred to Burgschweiger and Tsotsas [36]. Regarding the underlying assumptions of this model, we

summarise them as follows: (i) plug flow in the bubble phase; (ii) the particle-free bubble phase and the suspension phase within the bed are modelled separately, (iii) mass and heat transfer between drying gas and bubbles is significant and included in the model; (iv) heat transfer between the bed wall, particles, suspension gas, environment and bubble gas is also included.

#### *2.2. Deriving the Operational Envelopes*

As described in Samsatli et al. [26] the aim of deriving the operational envelopes of a process or unit operation is to find the maximum range of uncertain operating policies over which the design can be guaranteed to meet specific targets. The union of the maximum range of the uncertainty operating policies is referred to as the "operational envelope". This is particularly important for continuous pharmaceutical manufacturing as a multistage process, since through the use of such decoupled envelopes for each unit operation it can be ensured that the product specifications can be met if we restrict ourselves within the operating limits denoted through these envelopes.

The geometry of these envelopes can be arbitrary. However, in this work we employ hyperrectangular geometry for the sake of computational simplicity. Mathematically, if we denote by *b* ∈ - *bmin*, *bmax* the vector of uncertain parameters and their respective limits, which can be inferred either by expert knowledge or based on past observations, we seek to maximise the following objective function:

$$z = \prod\_{i=1}^{N\_b} b\_i^{\max} - b\_i^{\min} \tag{1}$$

where the index *i* = 1, ... , *Nb* is the index of the parameters under investigation. Instead of this objective function, which is non-convex, Samsatli et al. [26] proposed the use of a linear counterpart by introducing the difference in the magnitude of the ranges, i.e., Δ*bi* = *bmax i* − *bmin i* ∀*i*. Following this step, Equation (1) is replaced by the linear Equation (2) which reflects the scaled perimeter of the envelope.

$$f = \frac{1}{N\_b} \sum\_{i=1}^{N\_b} \frac{\Delta b\_i - \Delta b\_i^{\min}}{\Delta b\_i^{\max} - \Delta b\_i^{\min}} \tag{2}$$

Intuitively, since Equation (2) reflects a scaled perimeter the objective function range is [0,1] with an value of 0 reflecting the minimal envelope possible, i.e., <sup>Δ</sup>*bi* = Δ*bmin i* ∀*i*, and the maximal envelope feasible is obtained at the value of 1 where <sup>Δ</sup>*bi* = Δ*bmax i* ∀*i*. With this modification the overall problem that maximises *f* is given by model **(M1)**.

$$\begin{aligned} \max\_{a,b^{\min},b^{\max}} f &= \frac{1}{\Delta\_b} \sum\_{i=1}^{\Delta b\_i} \frac{\Delta b\_i - \Delta b\_i^{\min}}{\Delta b\_i^{\max} - \Delta b\_i^{\min}} \\ \text{Subject to} \\ \Phi\_0 \left[ \dot{\mathbf{x}}\_0, \mathbf{x}\_0, y\_0, a\_0, b\_0 \right] &= 0 \quad \forall b \in \left[ b^{\min}, b^{\max} \right] \\ h(\dot{\mathbf{x}}, \mathbf{x}, y, a, b) &= 0 \quad \forall b \in \left[ b^{\min}, b^{\max} \right], \ t \in (0, \tau], \ \tau \in \mathbf{b} \\ g(\dot{\mathbf{x}}, \mathbf{x}, y, a, b) &\ge 0 \; \forall b \in \left[ b^{\min}, b^{\max} \right], \ t \in (0, \tau], \ \tau \in \mathbf{b} \\ \Delta b &= b^{\max} - b^{\min} \\ \Delta b^{\min} &\le \Delta b \le \Delta b^{\max} \end{aligned} \tag{\mathbf{M1}}$$

In model **(M1)**, **Φ**0 represents the set of initial conditions for the system under study; *h*(·) represents the vector of equality constraints which are part of the model, e.g., mass/energy balances; *g*(·) represents the vector of inequality constraints, e.g., product specifications/resource limitations; *x* corresponds to differential state variables; *. x* their derivatives with respect to time (t); *y* represents algebraic state variables; while *a*, *b* represent time variant and time invariant controls, respectively. Notice that in **(M1)** the upper

bound of the time horizon is also allowed to be an "envelope" variable in case one wanted to investigate suitable bounds, for example for drying times.

Model **(M1)** is a semi-infinite programming problem since it needs to be solved for all the possible values of the *b* vector of variables. To overcome this issue, a twostep multiscenario optimisation problem is solved in which the envelope variables are discretised as described in Samsatli et al. [26].

#### **3. Case Study: Segmented Fluidised Bed Dryer**

In this section we demonstrate the methodology using the digital model of the continuous pharmaceutical process of the Diamond Pilot Plant (DiPP) at the University of Sheffield, shown in Figure 1. The process is a tableting pilot plant at the heart of which is a fluidised bed dryer (FBD) which is critical to the production of consistent quality product. The fluidised bed dryer (FBD) fluidises the feed granules to reduce their moisture content. In the process high-pressure hot air is introduced through a perforated bed of moist solid granules. The wet solids are lifted from the bottom and when fluidised are suspended in a stream of air. Heat transfer is accomplished by direct contact between the wet solid and hot gases. The vaporised liquid is carried away by the gas stream. The temperature and rate of input gas can be adjusted to save energy by, for example, aiming to shorten the drying time and manipulate the desired product (pharmaceutical granules) quality subject to a required range for the moisture content. The FBD is typically divided into a number of vertical segments.

As the FBD is connected with continuous twin screw granulation, the segmented FBD will ensure the wet granules in one cell are dried whilst the incoming wet granules flow into the neighbouring cell. Once the drying process in one cell is finished, the respective cell is emptied pneumatically and then conveyed to the downstream unit, in this case a mill. More segments contribute to reducing moisture but consume more time. In this study we set the FBD equipment to have two segments. Each segmen<sup>t</sup> size is 0.035 m3, with initial charge of 0.1 kg wet air and 0.1 kg granulates (lactose), with a particle density of 750 kg/m2. With these equipment specifications and initial conditions, the drying time is fixed by setting the volume and mass of the FBD, while temperature and flowrate of input streams are time-varied operating variables for achieving the moisture content objective. We implemented a single-factor experiment using gPROMS to investigate the effect of drying times and the two operational parameters, temperature and flowrate of input gas, on the envelope size. Using these studies enables us to find a suitable design that consumes less time and energy but has a bigger operational envelope.

Within a time interval *τ*0, *τf* , solid particles flow through cells of the FBD, and air with a temperature of *<sup>T</sup>*(*τ*) and a rate of *<sup>V</sup>*(τ) is continuously fed to the bottom of the FBD. Through fluidisation of the particles and consequent drying of the particles, the moisture content <sup>Γ</sup>(τ) of feed granules is reduced to the goal of a moisture content Γ (which could be a point or an interval). *V* is the volumetric flowrate and *T* is the temperature.

Employing the approach for traditional optimal control, we used the FBD model developed within gPROMS as a black box model [31], adding end point and path constraints. We used a black box model in order to show how it could be done without access to the full model equation set since this often needs to be the case in commercial settings.

The mathematical formulation is as follows:

> Γ(*t*) =

$$\min\_{\mathbf{x}, \mathbf{y}, T, V} f = \Gamma\_{\theta}$$

$$\text{Subject to : }$$

$$\Phi(\mathbf{x}(t), \mathbf{y}, T(t), V(t), \tau), 0 \le t$$

≤ *τ*

> *τ*]

(**M2**)

$$\text{with}$$

End point constraints : Γmin ≤ Γ*τ* ≤ ΓmaxΓ*τ* ≤ Γmax

$$\text{Path constraints : } \qquad T^{\text{min}} \le T(t) \le T^{\text{max}}, \forall t \in [0, 1]$$

*V*min ≤ *V*(*t*) ≤ *V*max, ∀*t* ∈ [0, *<sup>τ</sup>*],

where min and max refer to the upper and lower bounds, respectively, for each operational variable that is controllable. *x* and *y* refer to other model parameters that are uncontrollable. The drying time *τf*is a design variable and is fixed.

For each fixed value of the drying time, we applied the methodology shown in Section 2 to find an optimal operating envelope. We were then able to explore the design sensitivity by varying the value of the drying time to find a suitable design that consumes less time and energy but has a bigger operational envelope. The selected design would be the one that consumes less energy and has more flexibility.

Using the methodology shown in Section 2, to obtain an optimal balance between design and operational variables, we let *b* = -*T*min, *<sup>T</sup>*max, *V*min, *<sup>V</sup>*max, and formulate the following problem to determine the optimal operating envelope:

$$\max\_{y, b^{\min}, b^{\max}} f \equiv \frac{1}{N\_b} \sum\_{i=1}^{N\_b} \frac{\Delta b\_i - \Delta b\_i^{\min}}{\Delta b\_i^{\max} - \Delta b\_i^{\min}}$$
 
$$\text{Subject to : }$$

$$\Gamma'(\tau) = f(\mathbf{x}(\tau), y, b\_i, \tau), \qquad \tau\_0 \le \tau \le \tau\_f$$

$$\Gamma^{\text{min}} \le \Gamma\_{\tau\_f} \le \Gamma^{\text{max}} \text{ or } \Gamma\_{\tau\_f} \le \Gamma^{\text{max}} \tag{\text{M3}}$$

$$y^{\text{min}} \le b\_i \le y^{\text{max}}$$

$$\Delta b\_i = b\_i^{\text{max}} - b\_i^{\text{min}}$$

$$\Delta b\_i^{\text{min}} \le \Delta b\_i \le \Delta b\_i^{\text{max}}$$

The process modeling tool gPROMS [29] was used to implement and solve the model to determine the optimal operating envelopes. The gPROMS modeling platform allows existing models of processes to be converted to the envelope form and optimise their dynamic operation. The solution steps are briefly illustrated as follows:

Step 1: fix the value of design variable *τ*, the upper and lower bounds Δ*T*, Δ*V* and Γ, specify the interested range *T*min, *T*max *V*min, *V*max of the bounded variables, and let

$$T^{\text{min}} \le T^{\text{min}} \le T \le T^{\text{max}} \le T^{\text{max}}$$

$$\overline{V}^{\text{min}} \le V^{\text{min}} \le V \le V^{\text{max}} \le \overline{V}^{\text{max}}$$

Step 2: generate *NS* scenarios, each with a different set of operational variables ( *T*, *V*). For scenario *k* = 1, ··· , *Ns*, the values are given by:

$$T^{[i]} = T^{\min} + p^{[k]} \left( T^{\max} - T^{\min} \right)$$

$$V^{[i]} = V^{\min} + p^{[k]} \left( V^{\max} - V^{\min} \right)$$

where *p*[*k*] are normalized positions. For example, an optimization using two scenarios (*NS* = 2), one corresponding to the bottom left and another to the top right of the feasible region, we specify:

$$p^{|1|} = (0,0,\ldots,0), \ p^{|2|} = (1,1,\ldots,1)$$

Step 3: Then we define the objective function, variables and constraints from the FBD model within gPROMS, and solve the optimization problem to obtain the best values of *T*min, *T*max and *V*min, *V*max .

The algorithms were run on a personal computer with four 3.50 GHz processors and 16.0GB RAM using the Windows 10 operating system. The model and the approach can be used to optimise the steady-state and/or the dynamic behaviour of a continuous or batch process; in this case the fluid bed dryer is continuous.

The sampling technique employed in this work was a grid-based quasi-Monte Carlo sampling by using Sobol' low discrepancy sequences [37]. They have been shown to provide good distribution coverage even for fairly small sampling points. The design space was partitioned into a number of square grids and then within each grid sampling points were generated to evaluate feasibility. The interested reader is referred to Kucherenko et al. [38] for an in-depth discussion on the subject. In brief, for a response variable *<sup>Y</sup>*(*<sup>X</sup>*1, *X*2,..., *Xk*) which is a function of a set of input variables *X*1, *X*2, ... *Xk* a unit hypercube can be defined over the *k*-dimensions. Combining unit hypercubes over a grid-partitioned design space with quasi-random sequences is the most uniform possible solution to secure coverage. This is due to the fact that quasi-random points are selected from a sequence whilst knowing the position of the previous points and thus filling gaps between them [38].

We constructed an independent FBD model **(M2)**, to minimise drying time and moisture content, respectively, subject to it being in the interval [10%, 40%]. Next, we took the following steps:

Step 1: Specify the range of the operating variables:

$$\left( \left[ T^{\text{min}}, T^{\text{max}} \right] = \left[ 20 \text{ } \text{°C}, 80 \text{ } \text{°C} \right], \left[ V^{\text{min}}, V^{\text{max}} \right] = \left[ 240 \text{ m}^3/\text{h}, 480 \text{ m}^3/\text{h} \right] \text{h}$$

Step 2: Determine the feasible operating range with a drying time of 900 s which specifies a range of outputs of interest and hence a range of inputs. We uniformly sampled 13 temperatures in the range [20, 80] °C and 25 flow rates in the range [240, 480] m3/h. Next, we simulated the FBD model to detect the feasible region (i.e., 13 × 25 = 325 points) that satisfies end point and path constraints. Finally, we found all feasible solutions where the moisture falls in the range [10%, 40%]. This is shown in Figure 2.

Step 3: Run the optimisation model **(M3)** with a drying time of 900 s to obtain the operating envelope for *T* and *V*.


**Figure 2.** Feasible design range for *T* and *V* at a drying time of 900 s.

**Figure 3.** Operational envelope for a drying time of 900 s: *f* = 0.77.

**Figure 4.** Operational envelope for a drying time of 900 s while maintaining the maximal distance to the feasible boundary: *f* = 0.325.

The final stage is to explore the trade-off between design and operational flexibility as measured by the envelope size. The FBD model indicates that the feasible design space varies with the drying time. Hence, we can select a best drying time by exploring the envelope size. To do this we used a scenario-based algorithm with 10 candidate drying times (600–1500 s) and allowed Δ*T* and Δ*V* to vary.

From the results shown in Figure 5, we found that the FBD process can obtain the maximal envelope size with 700 s (as shown Figure 6 where a larger number of sampling points, i.e., 1000, was used to increase the resolution of the results), which means that this design has the best flexibility using the chosen operating variables. Figure 5 shows that there is significant effect on the flexibility of the process at different drying times with the optimal obtained at 700 s. Interestingly, in this case, the flexibility is not affected by the change in Δ*V* but only by the change in temperature, for the specified ranges of uncertainty. Nonetheless, we should point out that in this work the related nonlinear programming models were solved with a local and not a global optimisation solver which could explain some of the irregularities shown in Figure 5 for design options and envelope sizes.

**Figure 5.** Result of a design selection by trade-off between envelope size and drying time.

**Figure 6.** Operational envelope for a drying time of 700 s.
