**1. Introduction**

Pharmaceutical manufacturing processes have traditionally employed the batch operation mode, in which fixed amounts of raw materials are run through different unit operations to obtain the final drug product. Quality attributes of the final drug product were originally tested at the end of each batch processing step, where quality control essentially followed a quality-by-testing approach (QbT) [1], e.g., mixing uniformity is tested at the conclusion of the blending process. Over the last few years several factors have driven a shift from batch towards continuous pharmaceutical manufacturing. These factors include a reduction in the development cost for new medicines, making it both desirable and feasible to produce smaller annual volumes of targeted solutions for smaller patient populations, as well as improving product quality, decreasing cycle time, and better

**Citation:** Huang, Y.-S.; Sheriff, M.Z.; Bachawala, S.; Gonzalez, M.; Nagy, Z.K.; Reklaitis, G.V. Evaluation of a Combined MHE-NMPC Approach to Handle Plant-Model Mismatch in a Rotary Tablet Press. *Processes* **2021**, *9*, 1612. https://doi.org/10.3390/ pr9091612

Academic Editor: Luis Puigjaner

Received: 16 August 2021 Accepted: 3 September 2021 Published: 8 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>2</sup>School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907, USA;sbachawa@purdue.edu(S.B.);marcial-gonzalez@purdue.edu(M.G.)

controlled processing, to name a few popular drivers [2]. An economic analysis provided by Schaber and co-workers [3] highlights that continuous operation is able to provide estimated overall savings that can range from 9 to 40%, depending on the drug loading and process chosen, when compared to traditional batch operation. Additionally, encouraged by the regulatory agencies to modernize pharmaceutical manufacturing processes, academia and industry have invested significant time and resources to study different aspects required to successfully shift from batch to continuous operation mode. These efforts were made possible through various collaborations and consortiums [4–6].

In 2012, Gernaey and co-workers [7] identified the design and implementation of continuous pharmaceutical processes as one of the many issues that remain unresolved. Advanced process understanding is critical to the implementation of continuous pharmaceutical manufacturing applications [8]. To address this requirement, a quality-by-design (QbD) approach was pursued over the last decade [9]. QbD is a multi-step procedure that involves: (i) definition of quality target product profiles (QTPPs) and critical quality attributes (CQAs), (ii) identification of critical material attributes (CMAs) and critical process parameters (CPPs), (iii) linking of the CMAs and CPPs with the CQAs, (iv) examination of the design space and required control strategies, (v) validation, scale-up, and production [10]. While QbT primarily focused on end-stage testing, QbD revolved around advanced product and process understanding for systematic design of the operating space using mechanistic models and design of experiments (DoE). However, more recently there has been a shift towards quality-by-control (QbC), wherein quantitative and predictive understanding is leveraged for active process control with robust process design and operation, enabling smart manufacturing [11].

A goal of the QbC approach is real-time process monitoring and management, wherein advanced process control strategies are utilized to handle disturbances and exceptional events [11]. Process analytical technology (PAT) methods play a crucial role in monitoring a variety of CQAs in order to accomplish this [8,9]. Monitoring and control of CQAs such as tensile strength and tablet porosity are critical as they are linked to dissolution profiles of the manufactured tablets, which are ultimately linked to patient safety and treatment efficacy [12–16]. Tablet tensile strength and dissolution profile are affected by various factors such as particle size, API concentration, and addition of lubricants and glidant [17,18]. Glidants are added to improve the flowability of the blend. However, glidants and lubricants are also known to impact other product parameters, such as bulk density, compactibility, and compressibility. An objective of this work is to incorporate the impact of the use and control of glidants while assuring that the critical properties, such as tensile strength of the manufactured tablets, are maintained at desirable levels. In the context of continuous manufacturing, when a glidant feeder is used, it is important to use calibrated mechanistic models to handle the variations of glidant concentration. Therefore, it is essential to explore a variety of different control strategies to address the changes in CQAs of a tablet that may arise when glidants are used. It is worth noting that even though these challenges and control strategies are also relevant to lubricants, the development of mechanistic models and relevant control strategies associated with lubricants is beyond the scope of this work.

The identification and handling of plant-model mismatch (PMM) is an important component of any real-time process monitoring and control approach, and it has been an area of interest for decades. PMM can arise in the continuous tablet manufacturing process for a variety of reasons, e.g., the feeder refill step is known to introduce disturbances that affect CMAs such as the bulk density [19–21]. Since this may result in a deviation in the CQAs, PMM needs to be monitored and algorithms to mitigate it need to be developed and implemented. In order to monitor PMM, Harris initially presented a minimum variance-based assessment criterion to assess the condition of the working control loop [22]. This approach has gained popularity, but it is limited to single-input-single-output (SISO) systems [23]. More recently, data-driven methods that examine autocovariance and solve an optimization problem formulated to address the mismatch estimates in MIMO systems have been

developed. These methods minimize the discrepancy between the autocovariance of the output and the actual autocovariance of the mean-centered output variable [24–26]. Partial correlation coefficient (PCC)-based methods to identify PMM have also received attention in the literature, where PCC is well-suited as it is able to handle cases with high correlation in the manipulated variables [27,28]. As model re-identification is a critical, and often a time-consuming step once PMM has been identified, hybrid machine learning approaches have been proposed in order to aid model selection [28]. While there is grea<sup>t</sup> depth in the literature associated with the identification of PMM, there is limited discussion on practical approaches that would be applicable to the continuous pharmaceutical manufacturing industry in terms of managemen<sup>t</sup> of the PMM [29,30]. CQAs and CMAs need to be tracked online during plant operation but they may be unmeasurable in practice through existing PAT sensing methods (e.g., the bulk density within a unit operation); therefore, alternative solutions are required. This work proposes novel state estimation methods to accurately track states and model parameters online and, hence, guide operating decisions.

Additionally, most work in the continuous tablet manufacturing domain utilizes linear model predictive control strategies, often resulting from the linearization of a nonlinear system, which may not be adequate for some strongly nonlinear processes [31–33]. A literature review of traditional MPC application for different unit operations in the continuous pharmaceutical manufacturing industry can be found in [34], including end-to-end in silico and implementation studies. Since there is limited implementation of nonlinear model predictive control strategies for the continuous pharmaceutical manufacturing industry, a main objective of this work is to develop and present a moving horizon estimation-based nonlinear model predictive control (MHE-NMPC) framework to serve the dual requirement of accurate estimation and effective control. Model predictive control strategies are also advantageous over traditional proportional-integral-derivative (PID)-based control strategies, as they are able to effectively handle constraints, loop interactions and nonsquare control systems that may be encountered in manufacturing of pharmaceutical solid dosage forms [35–38]. A main practical concern for any developed framework is the need to ensure that the optimization problem can be solved in real-time, particularly for relatively quick processes such as those in the continuous pharmaceutical manufacturing industry. Therefore, an additional objective is to examine and discuss the real-time feasibility of the developed framework in controlling a rotary tablet press.

It is important to note that once non-conforming quality attributes have been identified, a long-term goal is the integration of control frameworks similar to the MHE-NMPC structure with residence time distribution (RTD)-based modeling frameworks that are currently being developed to guide tablet product diversion in the continuous pharmaceutical manufacturing industry and truly enhance and enable smart manufacturing operations [39,40].

To summarize, the primary objective of this work is to develop and present a moving horizon estimation-based nonlinear model predictive control (MHE-NMPC) framework to serve the dual requirement of accurate estimation and effective control, and to demonstrate its practical applicability by discussing its implementation feasibility in controlling a rotary tablet press. A secondary objective of this work is to examine different control strategies that are required when incorporating glidant feeders to further control tablet properties.

The rest of this work is organized as follows. In Section 2, mathematical modeling and optimization approaches for state estimation and control will be briefly discussed, along with the proposed monitoring algorithm and its advantages, i.e., the MHE-NMPC framework to monitor CMAs and CQAs and determine control actions will be presented. Section 3 will illustrate the robustness of the proposed MHE-NMPC framework with two examples of application. Specifically, the studies will showcase monitoring and control of a rotary tablet press in the presence of (i) plant-model mismatch and (ii) uncertainty in the glidant concentration. Section 4 will provide concluding remarks and directions for future work.

#### **2. Material and Methods**

*2.1. State Estimation*

State and parameter estimation methods have been utilized to enhance process monitoring capabilities in a number of industrial applications, ranging from bioreactors to robotics to continuous pharmaceutical manufacturing [30,41,42]. State and parameter estimation is a powerful tool in scenarios where process states or model parameters cannot be directly measured with sensors.

A nonlinear state-space model is defined as follows:

$$\mathbf{x} = \mathbf{g}(\mathbf{x}, \boldsymbol{u}, \boldsymbol{\theta}, \boldsymbol{w}) \tag{1}$$

$$y = l(\mathfrak{x}, \mathfrak{u}, \theta, \mathfrak{v}) \tag{2}$$

where *x*, *u*, *θ*, and *y* are the state variable vector, input variable vector, model parameter vector, and measurement vector, respectively [43]. The process and measurement noise are denoted by *w* and *v*, respectively. A schematic illustration of conventional state estimation algorithms is presented in Figure 1, where the nonlinear model is initialized based on the state values at the previous time step (*k* − 1) in order to obtain a prediction of the states and model parameters at the current time (*k*). State measurements are obtained from available sensors and are utilized to obtain a more accurate estimate of the states and parameter values by correcting the predictions from the model.

**Figure 1.** Conventional state estimation algorithm.

A number of alternative algorithms to carry out state estimation have been developed, e.g., the Kalman filter (KF), extended Kalman filter (EKF), unscented Kalman filter (UKF), particle filter (PF), and moving horizon estimation (MHE) are among the more popular approaches [44–47]. While the KF and EKF algorithms are suitable for linear or approximately linear applications, the UKF, PF, and MHE algorithms are able to handle processes that are more nonlinear in nature. The KF, EKF, and UKF algorithms also assume that error distributions are Gaussian in nature, while this assumption does not have to be satisfied for the PF and MHE algorithms [46]. Unlike the conventional approach illustrated in Figure 1, MHE utilizes a window, or moving horizon, of previous measurements in order to estimate the current states and model parameters, often providing improved performance when compared to the other algorithms. MHE is also capable of handling measurements collected from sensors at different sampling intervals or frequencies, which is advantageous for industries that utilize a variety of sensors to track physical attributes, as is the case in continuous pharmaceutical manufacturing. Therefore, MHE will be the choice of state estimation algorithm utilized in this work, as it is able to effectively track plant-model mismatch caused by deviations in model parameters, such as variations in the bulk density due to uncertainty in upstream unit operations (e.g., refilling of feeders).

The importance of monitoring powder feeder dosing in continuous pharmaceutical manufacturing is investigated in grea<sup>t</sup> detail by Destro and co-workers [29], wherein an MHE-based state estimation approach is implemented to reconcile mass measurements that are available from loss-in-weight (LIW) feeders with downstream measurements that are available from a PAT instrument and, thereby, to obtain practically continuous measurements as opposed to sampled measurements provided from the PAT instrument. The authors also demonstrate that the MHE approach is superior to one that utilizes statistical filters instead of the state estimator. Similarly, robust estimators were incorporated within the MHE skeleton of a feeding-blending system to handle dynamic systems with gross errors [30]. While the results presented in both case studies are promising, neither discusses or elaborates on the importance of their integration with efficient control strategies. To the knowledge of the authors, there has been no work in the continuous pharmaceutical manufacturing domain that has examined the integration of state estimation strategies with efficient control strategies, and therefore, this will be the central focus of this work, as discussed in the following sections.

#### *2.2. Model Predictive Control (MPC)—Linear and Nonlinear*

Model predictive control methods have been employed by various industries over the past few decades [48–52]. MPC relies on the dynamic model of the process. This model can either be linear or linearized models obtained through system identification, as in the case of the linear implementation of MPC, or be nonlinear and derived using first principles or semi-empirically (using a hybrid model), as in the case of NMPC [53]. Both MPC algorithms utilize a finite time horizon to optimize the control input at the current time iteration, while keeping future time iterations in mind. This ability makes MPC predictive in nature due to its ability to anticipate future events and take control actions accordingly, which is not possible using traditional PID controllers [54].

While nonlinear model predictive control (NMPC) methods have been utilized by some industries, to the knowledge of the authors their implementation has not been explored extensively for continuous pharmaceutical manufacturing of solid dosage forms, cf. [31–33], where linear or hybrid implementations of MPC are utilized. Since processes in the continuous pharmaceutical manufacturing industry are known to be nonlinear in nature, it is therefore desirable to develop and implement an NMPC approach.

It should also be noted that these predictive control strategies are particularly advantageous for cases of non-square systems, i.e., where the number of manipulated variables exceeds the number of the controlled variables, since these methods are able to effectively manage nonlinear relationships [35–37]. These cases cannot be straightforwardly handled using traditional PID control strategies [38].

The following section will present the developed MHE-NMPC framework that seeks to accomplish the dual requirement of accurate estimation and efficient control. Since real-time implementation feasibility is an area of interest, a discussion on the practical applicability of the framework developed will also be presented.

#### *2.3. Moving Horizon Estimation-Based Nonlinear Model Predictive Control (MHE-NMPC) Framework*

The algorithm proposed in this work seeks to combine the effective estimation capabilities of MHE with the control abilities of NMPC, through the MHE-NMPC framework illustrated in Figure 2. Real-time measurements of output variables (*y*) and input variables (*u*) are first collected to monitor the process. Since disturbances, either known or unknown, can always exist in a real plant, mismatches may arise between the sensor measurements and model values. As elaborated previously, the goal of state estimation is to obtain a 'true state' value by utilizing the information from both measurements and process models. The 'true state' can be either measurable, e.g., API concentration at the blender exit using NIR sensors [55–57], or unmeasured, e.g., powder holdup in the blender. Through the updating of uncertain model parameters, which have changes due to upstream disturbances, MHE enables the handling of plant-model mismatch, thus allowing the controller to receive estimated output variables (*y*<sup>ˆ</sup>) with less uncertainty. The NMPC control algorithm then minimizes the error between setpoints *ysp* and estimated output variables (*y*<sup>ˆ</sup>) by deciding the optimal control move (*u*) for the process to reach both setpoint tracking and

disturbance rejection, i.e., the control objectives, while updating the model parameter (<sup>ˆ</sup>*θk*) and median of the error distribution in the past time window (*ζ*).

**Figure 2.** Adaptive control framework of MHE-NMPC.

A schematic illustration of the MHE-NMPC framework at time *t* = *k* is shown in Figure 3, with *Npast* measurements available in the past window and *Np* estimations in the prediction window. The MHE is then formulated as follows:

$$\min\_{\hat{\theta}\_k} J = \sum\_{t=k-N\_{\text{post}}}^k (\mathfrak{e}\_t)^T \mathcal{W}\_E \,\mathfrak{e}\_t + (\hat{\theta}\_k - \hat{\theta}\_{k-1})^T \mathcal{W}\_\theta (\hat{\theta}\_k - \hat{\theta}\_{k-1}) \tag{3a}$$

subject to

$$\pounds\_{k-N\_{pust}+j+1} = f(\pounds\_{k-N\_{pust}+j\_{\prime}}\mu\_{k-N\_{pust}+j\_{\prime}}\,\hat{\theta}\_{k})\tag{3b}$$

$$\mathfrak{F}\_{k-N\_{pust}+j} = h(\mathfrak{x}\_{k-N\_{pust}+j}) \tag{3c}$$

$$
\varepsilon\_{k-N\_{\text{post}}+j} = \mathcal{Y}\_{k-N\_{\text{post}}+j} - \mathcal{Y}\_{k-N\_{\text{post}}+j} \tag{3d}
$$

$$\mathfrak{A}\_{k-N\_{\rm post}+j+1} \in \mathbb{X}\_{\prime} \quad \mathfrak{e}\_{k-N\_{\rm post}+j} \in \Omega\_{\mathbb{C}} \quad \hat{\theta}\_{k} \in \Omega\_{\theta} \tag{3e}$$

$$j = 0, \ 1, \ \ldots, \ N\_{post} \tag{3f}$$

where ˆ *θk* are estimated uncertain parameters, which are bounded in the compact set Ω*θ*. In the above formulation, *yt* and *ut* are measurements of output variables and input variables at time *t*, respectively; *y*ˆ*t* and *x*ˆ*t* are estimated output and state values, respectively; *εt* are output disturbances, which are bounded in the compact set Ω; and *WE* and *Wθ* are the weighting matrices. After the MHE optimization problem is solved at time *t=k*, the estimated state *<sup>x</sup>*<sup>ˆ</sup>*k*−*Npast*+<sup>1</sup>|*<sup>t</sup>*=*<sup>k</sup>* is chosen as the initial state value of next time step *t=k+* 1, i.e., *<sup>x</sup>*<sup>ˆ</sup>*k*−*Npast*+<sup>1</sup>|*<sup>t</sup>*=*k*+<sup>1</sup> = *<sup>x</sup>*<sup>ˆ</sup>*k*−*Npast*+<sup>1</sup>|*<sup>t</sup>*=*<sup>k</sup>* [58].

While an error distribution of output variables *yt* − *y*ˆ*t* in the past time window can be obtained from Equation (3d), a single point estimate of the output *y*ˆ*t* is of most interest in many applications, instead of the whole error distribution [59]. When no probabilistic process models are used, it is easier to use a single point estimate of the output *y*ˆ*t* to visualize and control process dynamics. In this study, the median of the error distribution in the past time window is used to represent output disturbances *ζk*at time *t* = *k*, i.e.,

$$\mathcal{Z}\_k = median \begin{Bmatrix} \varepsilon\_{k - N\_{\text{post}} + j} \end{Bmatrix}, \text{for } j = 0, \ 1, \dots, N\_{\text{post}} \tag{4}$$

Therefore, with estimated states *<sup>x</sup>*<sup>ˆ</sup>*k*, output disturbances *ζ<sup>k</sup>*, and updated uncertain optimal parameters ˆ *θk*, the NMPC framework at time *t* = *k* is given by:

$$\min\_{\Delta u\_t} J = \sum\_{t=k}^{k+N\_p} \left( \mathcal{y}\_t - y\_{sp} \right)^T \mathcal{W}\_y \left( \mathcal{y}\_t - y\_{sp} \right) + \sum\_{t=k}^{k+N\_c-1} \left( \Delta u\_t^T \mathcal{W}\_{\Delta u} \Delta u\_t \right) \tag{5a}$$

ˆ

subject to

$$\mathfrak{k}\_{k+j+1} = f(\mathfrak{k}\_{k+j\prime} \hat{u}\_{k+j\prime} \, \theta\_k) \tag{5b}$$

$$\mathfrak{F}\_{k+j} = h \begin{pmatrix} \mathfrak{x}\_{k+j} \\ \end{pmatrix} + \mathfrak{J}\_k \tag{5c}$$

$$
\Delta u\_{k+j} = u\_{k+j+1} - u\_{k+j} \tag{5d}
$$

$$\Lambda \pounds\_{k+j} \in \mathbb{X}\_{\prime} \quad u\_{k+j} \in \mathbb{U}\_{\prime} \quad \Delta u\_{k+j} \in \Omega\_{\Lambda u} \tag{5e}$$

ˆ

$$j = 0, \ 1, \ \ldots, \ N\_p - 1 \tag{5}$$

ˆ

ˆ

ˆ

where *Nc* is the length of the control time window, and *ysp* are the setpoints of the output variables. *Wy* and *W*Δ*u* are the weighting matrices. Control movements Δ*u* are constrained in the compact set ΩΔ*<sup>u</sup>*. The control window *Nc* is usually smaller than the prediction window *Np* and has to be chosen considering a compromise between computational burden and stability requirements. Control movements <sup>Δ</sup>*uk*+*<sup>j</sup>* in control window *Nc* vary according to results of optimization, but those beyond the control window are zero, i.e., <sup>Δ</sup>*uk*+*Nc*=

ˆ

ˆ

$$
\Delta u\_{k+N\_\varepsilon+1} = \dots = \Delta u\_{k+N\_p-1} = 0,\\
\text{which implies that } u\_{k+N\_\varepsilon} = u\_{k+N\_\varepsilon+1} = \dots = u\_{k+N\_p}. \text{ In } \varepsilon \text{ we have } \varepsilon \ge 0, \text{ and } \varepsilon > 0, \text{ and } \varepsilon > 0.
$$

other words, while the predicted *<sup>y</sup>*<sup>ˆ</sup>*k*+*<sup>j</sup>* can still be calculated using <sup>Δ</sup>*uk*+*<sup>j</sup>* and *uk*+*j* in the prediction window *Np*, only <sup>Δ</sup>*uk*+*<sup>j</sup>* in control window *Nc* is considered in the objective function. It should be noted that models of estimated output variables *y*ˆ*t* are different in Equations (3c) and (5c). In the future time window, output disturbances *ζk* are added to the model of the process, allowing zero steady-state offset in controlled output variables *y* [59,60]. A schematic illustration of MHE-NMPC is provided in Figure 3, where at each iteration, the MHE is utilized to obtain a more accurate representation of the true state of the process and plant-model mismatch, and the NMPC is utilized to find the optimal first move for each input variable *u*. This framework thus allows for both accurate estimation and efficient control.

**Figure 3.** Illustration of MHE-NMPC coupling at each time interval.

#### *2.4. Implementation of a Real-Time Feasible MHE-NMPC Framework*

The MHE-NMPC framework is implemented in MATLAB (MathWorks R2018a) and the MATLAB built-in *fmincon* function is used to solve the optimization problem in each iteration. The computation is performed on a 64-bit ASUS VivoBook with Intel® Core ™ i7-8550U @1.80 GHz processor and 8GB of total memory. In all simulated results, the time unit for each step is 1 s, the past time window *Npast* used in MHE is 30 time units, and the NMPC is tuned with prediction time window *Np* chosen to be 60 time units and control time window (*Nc*) to be 10 time units. Sensor measurements are also assumed to be available at 1 s intervals. It should be noted that the average computation time for each iteration is 0.7 s, indicating that the optimization problem can be solved and implemented in real time. These results demonstrate the feasibility of the proposed framework, and its ability to achieve real-time process control.

The following section will explore the applicability of the developed MHE-NMPC framework to track plant-model mismatch and to efficiently control a key process unit operation in the continuous manufacturing line of solid dosage forms, i.e., the rotary tablet press.

#### **3. Examples of Application to Continuous Direct Compression**

The applicability of the developed MHE-NMPC framework will be demonstrated through two case studies. The first case study will highlight the importance of monitoring model parameters in real time and how this is enabled via state estimation, as opposed to the use of fixed model parameters. Different degrees of plant-model mismatch will be used. The second case study will present the applicability of the framework in the practical scenario of having uncertainty in the glidant feeding conditions.

Both case studies will focus primarily on the tablet press unit operation of the direct compression line. A hierarchical implementation of the three-level quality-by-control (QbC) framework of control systems for the continuous direct compression line is illustrated in Figure 4, whose unit operations are comprised of feeders, blenders, and the tablet press. For this line, Level 0 control is generally implemented through programmable logic control (PLC) systems built into the unit operation equipment in order to control CPPs. Level 1 control relies on PAT tools to monitor and control CQAs and can encompass multiple unit operations designed to reduce the impact of disturbances that may propagate downstream. Level 1 control supervises the Level 0 control, typically accomplished through SISO control loops which aim to maintain desired setpoints for the CQAs. A distributed control system (DCS) is employed to integrate process equipment such as the feeders and tablet press and any instrumentation that measures material properties. More advanced approaches are applied at Level 2 and use mathematical models such as MHE for validating process measurements, with the ability to predict the effects of disturbances and changes in the CPPs on the CQAs. Additional functionalities provided at Level 2 can include NMPC, a quality managemen<sup>t</sup> system (QMS), and real-time optimization (RTO).

**Figure 4.** A 3-level hierarchical implementation of control systems for the continuous direct compression process (modified from source: [11]).

#### *3.1. Tablet Press Model*

The rotary tablet press and the lubricant/glidant feeder are key unit operations, where the latter is used to reduce frictional loses and facilitate powder flow during die filling and formation of solid tablets via mechanical compression. Therefore, models for glidant effects in die filling and compression processes will be used to monitor and control the porosity and tensile strength of tablets. Specifically, these mechanistic models capture the effects of glidant concentration and mixing conditions [61,62].

The weight of a convex tablet, *W*, formed using Natoli D-type tooling with shallow cup depth can be computed as follows:

$$\mathcal{W} = \rho\_b V\_{fill} \left( 1 - \zeta\_1 \frac{n\_T}{n\_F} + \zeta\_2 \frac{H\_{fill}}{D} \right) \tag{6}$$

where *D*, *Vfill*, *Hfill*, *ρb*, *nT*, and *nF*, refer to the diameter of the die, volume of the die cavity, dosing position, bulk density of the powder, turret speed, and feed frame speed, respectively [62]. In Equation (6), *ξ*1 and *ξ*2 are model parameters to be estimated from experimental data. The bulk density depends on glidant concentration and mixing conditions, but its dependency is beyond the scope of this work. For the D-type tooling, the volume of the die cavity is given by:

$$V\_{fill} = \frac{\pi D^2 H\_{fill}}{4} + \frac{\pi h \left(\frac{3D^2}{4} + h^2\right)}{6} \tag{7}$$

where *h* is the cup depth. The tablet production rate, .*mtablet*, is given by:

$$
\dot{m}\_{tablet} = \mathcal{W}n\_T\mathcal{N}\_{tation} \tag{8}
$$

where *Nstation* is the number of turret stations available. For a blend composed of MCC (Avicel PH200), APAP (acetaminophen) and silica experimental evidence suggests that both pre-compression and main compression forces do not show a dependency on glidant conditions [62]. The pre-compression force (PCF) is then computed as follows [63]:

$$F\_{\rm pc} = \frac{\pi D^2}{4b} \left[ \frac{\rho^{\rm pc} - \rho\_c}{\rho^{\rm pc} (a - 1) + \rho\_c} \right] \tag{9}$$

where parameters *a* and *b* are Kawakita constants [63], which represent the maximum degree of compression and the reciprocal of the pressure applied to attain this degree of compression, respectively. In Equation (9), *ρc* and *ρpc* refer to the critical density and the pre-compression relative density, respectively. The pre-compression relative density is computed as follows:

$$
\rho^{\rm pc} = \frac{\mathcal{W}}{V^{\rm pc} \rho\_t} \tag{10}
$$

and

$$V^{pc} = \frac{\pi D^2 H^{pc}}{4} + \frac{\pi h \left(\frac{3D^2}{4} + h^2\right)}{3} \tag{11}$$

where *ρt* and *Hpc* refer to the true density of the powder and the pre-compression thickness, respectively. Similarly, the main compression force (*Fpunch*) is computed as follows:

$$F\_{\rm punch} = \frac{\pi D^2}{4b} \left[ \frac{\rho^{in - di\varepsilon} - \rho\_c}{\rho^{in - di\varepsilon} (a - 1) + \rho\_c} \right] \tag{12}$$

with the in-die relative density *ρin*−*die* given by:

$$
\rho^{\dot{m}-\dot{d}\dot{c}} = \frac{\mathcal{W}}{V^{\dot{m}-\dot{d}\dot{c}}\rho\_t} \tag{13}
$$

and

$$V^{in-\text{dir}} = \frac{\pi D^2 H^{in-\text{dir}}}{4} + \frac{\pi h \left(\frac{3D^2}{4} + h^2\right)}{3} \tag{14}$$

where *Hin*−*die* refers to the main compression thickness. The tablet density, or out-of-die relative density of the tablet, *<sup>ρ</sup>tablet*, is then obtained utilizing the elastic recovery, *ερ*, of the tablet as follows:

$$
\rho^{\text{tablet}} = (1 - \varepsilon\_{\rho}) \rho^{in - \text{dis}} \tag{15}
$$

The elastic recovery model is not sensitive to the glidant mixing conditions [61,62], and it is governed by:

$$
\varepsilon\_{\rho} = \varepsilon\_0 \frac{\rho^{\rm in-dir} - \rho\_{\rm c,c}}{1 - \rho\_{\rm c,c}} \tag{16}
$$

where *ε*0 and *ρ<sup>c</sup>*,*<sup>ε</sup>* are the in-die elastic recovery at full compaction and the relative density at which tablets do not exhibit elastic recovery, respectively [64]. The tensile strength *σt* exhibits dependency on glidant conditions and it is computed as follows

$$
\sigma\_t = \sigma\_0 \left[ 1 - \left( \frac{1 - \rho^{tabllet}}{1 - \rho\_{\varepsilon, \sigma\_t}} \right) e^{(\rho^{tabllet} - \rho\_{\varepsilon, \sigma\_t})} \right] \tag{17}
$$

where *σ*0 is the tensile strength at zero porosity and *ρ<sup>c</sup>*,*σ<sup>t</sup>* is the critical relative density at which tablets do not exhibit any the tensile strength, i.e., the relative density at which a tablet starts forming [17]. It bears emphasis that these parameters are functions of glidant conditions, specifically:

$$
\rho\_{\mathfrak{c},\mathcal{O}\_{\mathbb{T}}} = \frac{\rho\_{\mathfrak{c},0} - \rho\_{\mathfrak{c},\infty}}{1 + \mathbb{C}\_{\mathfrak{c}}} + \rho\_{\mathfrak{c},\infty} \tag{18}
$$

$$
\sigma\_0 = \frac{\sigma\_{0,\phi}}{1 + \mathbb{C}\_{\sigma}} \tag{19}
$$

$$\mathcal{C}\_{\sigma} = \frac{c\_I^{b\_1} \gamma^{b\_2}}{b\_3} \tag{20}$$

where *ρ<sup>c</sup>*,0, *ρ<sup>c</sup>*,∞, *<sup>σ</sup>*0,*φ*, *b*1, *b*2, and *b*3 are model parameters estimated from experimental data. In Equation (20), *cl* and *γ* are the glidant concentration and total shear imparted to the blend, respectively. For simplicity, the total shear strain is represented by an equivalent mixing time, which, in turn, is estimated as follows

$$
\gamma = \gamma\_0 + \frac{m\_{f,h}}{\dot{m}\_{\text{tablelet}}} \tag{21}
$$

where *γ*0 and *mf* ,*h* are a total shear strain base line, expressed in term of mixing time, and the mass hold up in the feed-frame and hopper, respectively. Specifically, for the purpose of these case studies, the model parameter *γ*0 is the glidant mixing time used in the rotary Tote blender when the blend was prepared. A 5L rotary Tote blender was employed. The mean residence time in the feed-frame, i.e., *mf* ,*h* .*mtablet* , is used to estimate the additional shear, or mixing time, imparted inside the tablet press.

The dependency of the bulk density on the glidant concentration, *cl*, can be incorporated through the following equation [62]:

$$
\rho\_{\mathbb{b}} = \rho\_{\mathbb{b},\infty} - \frac{\rho\_{\mathbb{b},\infty} - \rho\_{\mathbb{b},0}}{1 + \mathbb{C}\_{\rho}} \tag{22}
$$

where *ρb*,<sup>∞</sup> and *ρb*,<sup>0</sup> represent the bulk densities when the shear strain imparted is infinite and zero, respectively. *Cρ* is a lumped parameter that defines the glidant mixing conditions computed as follows [62]:

$$\mathbf{C}\_{\rho} = \frac{c\_1^{r\_1} (\gamma + \gamma\_0)^{r\_2}}{r\_3} \tag{23}$$

where *γ* and *γ*0 are the shear imparted to the powder during mixing and the initial shear strain imparted prior to mixing, respectively. *r*1, *r*2, and *r*3 are fitting parameters.

A Natoli NP-400 tablet press and SOTAX AT4 tablet tester were used in this work to fabricate tablets and gather experimental data under steady-state conditions. The experimental data were then used to carry out parameter fitting using the *fmincon* function in MATLAB to obtain realistic model parameters values that could be used for the simulations presented in case studies 1 and 2, which are summarized in Table 1.




**Table 1.** *Cont.*

*3.2. Case Study 1: Monitoring and Control of the Rotary Tablet Press in the Presence of Plant-Model Mismatch*

Monitoring powder bulk density in the tablet press is of critical importance, as it affects the tablet properties [12]. Sources of variability can be introduced during any of the unit operations upstream, e.g., in the feeder unit operations during refill, as the feeder switches from gravimetric mode to volumetric mode, leading to either increases in bulk density due to compression or decreases in bulk density due to aeration [19–21].

For this case study, a four-by-five system was employed as it would enable the incorporation of an extra manipulated input for added control benefits, i.e., glidant concentration. It is assumed that the direct compression line has the ability to utilize the glidant concentration as a manipulated variable through changes in the glidant flowrate. In practice this would be implemented in the hierarchical three-level QbC framework, by using a level-one PID control, that would use the glidant concentration measurement and adjust the glidant flowrate to follow the concentration setpoint set by the level-two NMPC. In this case the four-by-five non-square level-e control system is comprised of controlled variables consisting of the tablet weight, pre-compression force, production rate, and tensile strength and manipulated variables consisting of the dosing position, pre-compression thickness, main compression thickness, turret speed, and glidant concentration. It is assumed that measurements for the tablet weight, pre-compression force, main compression force, and production rate are all available every second [61]. In this simulation, it should be noted that the main compression force was not a directly controlled variable with specified set points. This is because the tensile strength could not be maintained while simultaneously fixing the main compression force. Given the objective to maintain the tensile strength at desired levels due to its link to patient safety, the tensile strength was chosen over the main compression force as a controlled variable. Since measurements of the main compression force are available, they were utilized in the MHE framework only for the purpose of parameter estimation. As maintaining the CQAs is important, and since the tensile strength

measurements are not available in real time, a soft sensor based on Equation (17) is utilized in order to track this particular state in real time. In practice, the SOTAX AT4 tablet tester can be utilized in order to obtain measurements of the tensile strength. However, since the diametrical compression test is destructive, tensile strength measurements are available at a lower frequency than one of the PAT sensors, which in turn also drives the need for a soft sensor.

A summary of the controlled variables, manipulated variables, measured variables, and uncertain model parameter is provided in Table 2. In order to examine the performance of the MHE-NMPC framework under PMM, three different scenarios will be examined, namely: nominal operation (no PMM), operation with mild PMM, and operation with high PMM. A summary of model parameters for each scenario is provided in Table 1, where the MHE-NMPC tuning parameters are those described in detail in Section 2.4. Mild PMM is simulated by introducing mismatch to three model parameters: *ρb*, *ρc* and Kawakita parameter *a*. High PMM is simulated by introducing mismatch to six model parameters: *ξ*2, *ρb*, *ρ<sup>c</sup>*, Kawakita parameters *a* and *b*, and *ρ<sup>t</sup>*. In this simulation, the 'model' and 'plant' share the same equations detailed in Section 3.1. Different parameter values were assigned to the 'model' and 'plant' in order to simulate mismatch. Additionally, sensor measurement noise in the plant was simulated by adding normally distributed error with zero mean and variance analogous to the variability of a real sensor. This variability was obtained from historical plant data.

**Table 2.** Summary of variables and uncertain model parameters for case study 1.


For all three scenarios, three model parameters are tracked, i.e., *ρb*, *ρc* and Kawakita parameter *a*. The bulk density was monitored due to its influence on a number of other model parameters and states, while the relative critical density and Kawakita parameter *a* are both known to influence the compression forces, making them critical parameters that also need to be tracked in real time.

Simulation results of the process outputs for all three scenarios are presented in Figure 5a–c, respectively. The MHE-NMPC framework is utilized for all simulations and includes open-loop control from 0–100 s (indicated by red shading in all plots), state estimation using MHE and open loop control from 100–200 s (indicated by yellow shading in all plots), and implementation of the MHE-NMPC framework from 200 s until the end of the simulation (indicated by gray shading in all plots). Setpoint changes are introduced for the tablet weight from 210 mg to 240 mg at 400 s, for the pre-compression force from 0.3 kN to 0.6 kN at 600 s, for the production rate from 6.9 kg/h to 8 kg/h at 800 s, and for the tensile strength from 4.2 MPa to 6 MPa at 600 s, respectively.

**Figure 5.** Controlled variables for case study 1 under scenarios with (**a**) no PMM, (**b**) mild PMM, and (**c**) high PMM.

For the scenario where there is no PMM, accurate setpoint tracking can be achieved for all states (see Figure 5a). This is also true for the case where there is mild PMM (see Figure 5b). This is primarily enabled due to the ability of the MHE-NMPC framework to accurately track variations in the uncertain model parameters in real time as illustrated in Figure 6b, allowing the impact of PMM to be effectively managed. In the case where there is high PMM (see Figure 5c), fairly accurate setpoint tracking can still be achieved despite there being mismatch in more parameters than those being tracked. This is also observed from the parameter estimation results in Figure 6c, although there is a slight offset in the model parameters being tracked to compensate for the variation in the additional model parameters that are not being tracked. The corresponding plots of the manipulated variables for all three scenarios are presented in Figure 7.

**Figure 6.** Parameter estimation of model parameters for case study 1 under scenario with (**a**) no PMM, (**b**) mild PMM, and (**c**) high PMM.

Control performance metrics are required in order to assess the performance of the development framework and accurately assess the impact of PMM for each scenario. Beyond the typical controller performance metric, i.e., integral of absolute error (IAE), some additional metrics are used to quantify the control performance in each scenario, which include the duration-to-reject (D2R) and magnitude to product (M2P) [65]. These metrics are able to quantify control performance in a manner that is more relevant for the continuous pharmaceutical industry. D2R is the duration of time that the process requires to smooth out the process disturbance or to reach a new set point for the CQA/CPP. M2P describes the maximum deviation in the CQA/CPP from the target setpoint. Larger values of all these performance metrics indicate worse or degraded control performance. The IAE values are calculated from t = 300 s to t = 1000 s. A summary of the control performance metrics is provided in Table 3.

When mild PMM exists, the control performance of the MHE-NMPC framework is comparable to the scenario without PMM, implying that the framework is able to sufficiently handle the PMM. However, when there is excessive PMM in multiple parameters as in the scenario with high PMM, significantly higher values of IAE and M2P are obtained, particularly in the tensile strength. This can be attributed to the fact that there was mismatch in more model parameters than those being tracked for this particular simulation, as can be noted from Table 1, resulting in an offset in the estimates of the model parameters to compensate for the added uncertainty (see Figure 6c). It should also be noted that the setpoints for the tablet weight and production rate track reasonably well, even in the

presence of high PMM as demonstrated by the comparable IAE values for both states. However, since the tensile strength is an important CQA that is linked to the dissolution profile of the tablets, once high PMM begins to cause an offset in the tensile strength, it can serve as an indicator for the requirement to carry out model re-identification. This case study was able to demonstrate the strength of the MHE-NMPC framework in its ability to handle PMM in multiple model parameters.

**Figure 7.** Manipulated variables for case study 1 under scenarios with (**a**) no PMM, (**b**) mild PMM, and (**c**) high PMM.


**Table 3.** Control performance of the MHE-NMPC framework for case study 1 under different levels of PMM.

#### *3.3. Case Study 2: Monitoring and Control of the Rotary Tablet Press in the Presence of Uncertainty in the Glidant Concentration*

This scenario aims to explore a more practical concern with regards to the incorporation of the glidant feeder in the control scheme. In practice in some applications, it might not be possible to accurately control the concentration of the glidant in the direct compression process, due to the low concentrations used and, thus, the small feeding rates needed. Uncertainty in the glidant concentration is important, as it leads to variations in the bulk density upstream of the rotary tablet press. Therefore, accurate monitoring and control of these variations is required.

Since the glidant concentration can no longer be treated as a manipulated input for this scenario, the original system is revised to form a four-by-four MIMO system with the controlled variables consisting of the tablet weight, pre-compression force, production rate, and tensile strength. The manipulated variables consist of the dosing position, precompression thickness, main compression thickness, and turret speed. Once again it is assumed that only measurements for the tablet weight, pre-compression force, main compression force, and production rate are available every second. As tensile strength measurements are not available in real time, a soft sensor based on Equation (17) is once again utilized for this particular state. The concentration of silica is then assumed to be an uncertain parameter. A summary of the controlled variables, manipulated variables, measured variables, and uncertain model parameters is provided in Table 4. A summary of the model parameters was provided in Table 1, where the MHE-NMPC tuning parameters are those detailed in Section 2.4.


**Table 4.** Summary of variables and uncertain model parameters for case study 2.

For this particular case study, mismatch is introduced through positive and negative step changes in the silica concentration from its nominal value of 0.2% to 0.35% between 300 and 700 s and from 0.2% to 0.05% between 1100 and 1500 s, respectively. Step changes in either direction are introduced in order to examine if and how the direction of the mismatch affects the control performance. Simulation results of the process outputs for open loop control, closed loop control using only NMPC control, and closed loop estimation and control using the proposed MHE-NMPC framework are presented in Figure 8a–c, respectively. The simulation for the NMPC framework includes open loop control from 0

to 200 s, with a closed loop NMPC strategy implemented from 200 s until the end of the simulation. The simulation for the MHE-NMPC framework includes open loop control from 0 to 100 s, state estimation using MHE and open loop control from 100 to 200 s, and implementation of the MHE-NMPC framework from 200 s until the end of the simulation.

**Figure 8.** Controlled variables for case study 2 under scenarios using (**a**) open loop control, (**b**) only NMPC, and (**c**) MHE-NMPC.

Figure 8a demonstrates failure of the open loop control to maintain the controlled variables at their setpoints and, thus, the need to implement effective control strategies. Due to the disturbance terms utilized in the NMPC framework, offset free control is achieved for three of the four controlled variables when only NMPC is employed, as illustrated in Figure 8b. However, since real-time measurements are unavailable for the tensile strength, a soft sensor is employed. Since the soft sensor does not incorporate a disturbance term, an offset can be observed between the soft sensor values and the set point due to the mismatch caused by the assumption of fixed model parameters in the NMPC framework. In contrast, as illustrated in Figure 8c the MHE-NMPC framework is able to achieve offset free control for all four states, including the tensile strength. This is attributed to the ability of the MHE-NMPC framework to provide (i) real-time and accurate estimation of uncertain model parameters, enabled by MHE, and (ii) efficient control, enabled by NMPC.

Figure 9 shows results for the uncertain model parameter estimation, i.e., the estimation of silica concentration. These results demonstrate the ability of the MHE-NMPC framework to accurately track variations in silica concentration. It should be noted that the sluggish behavior at higher concentrations of silica is due to the nonlinear effect silica has on the process. Since in practice the true value of the concentration of silica might be unknown, this case study demonstrates the advantage of utilizing the MHE-NMPC framework to achieve accurate estimation of both measurable and unmeasurable states and parameters, such as the concentration of silica, while executing realistic and effective control strategies.

**Figure 9.** Parameter estimation for case study 2 under scenarios using (**a**) open loop control, (**b**) only NMPC, and (**c**) MHE-NMPC.

Figure 10 shows the manipulated variables for the different control strategies studied in this section. All changes in the manipulated variables presented in these results are realistic in nature and can be achieved during normal operation of the tablet press. It should be noted that the variations in the manipulated variables are larger in the case where only NMPC is utilized (see Figure 10b) compared to the case where the MHE-NMPC framework is utilized (see Figure 10c). This may be attributed to a less effective linear output disturbance model implemented in the NMPC framework, when compared to the MHE-NMPC framework that also carries out parameter updating incorporating more directly nonlinear effects of the PMM in the scheme.

This case study was able to demonstrate the ability of the MHE-NMPC framework to track and manage fluctuations in the glidant concentration, often caused by upstream disturbances, thereby providing an efficient solution to a common process upset faced when operating the rotary tablet press.

**Figure 10.** Manipulated variables for case study 2 under scenarios using (**a**) open loop control, (**b**) only NMPC, and (**c**) MHE-NMPC.

### **4. Conclusions**

The continuous pharmaceutical manufacturing industry is in need of improved realtime process monitoring and managemen<sup>t</sup> strategies that, specifically, are able to effectively identify and handle plant-model mismatch (PMM). In order to enable the quality-bycontrol (QbC) paradigm to move forward, this work developed and presented a moving horizon estimation-based nonlinear model predictive control (MHE-NMPC) framework to accomplish the dual requirement of accurate estimation and efficient control. Realtime implementation feasibility of the developed framework was also discussed, and the ability of the proposed framework to solve the optimization problem at each time step in a manner that enabled real-time implementation was highlighted. The practical applicability of the developed framework was corroborated through two realistic case studies that incorporated the effects of glidant to better control CQAs such as the tensile strength. Both examples demonstrated the ability of the framework to achieve reasonable control performance despite the presence of varying sources and degrees of PMM.

Future work includes further demonstration of the practical applicability of the proposed MHE-NMPC framework utilizing the rotary tablet press at Purdue University, including the application of the framework to the entire direct compression line. While a soft sensor was utilized in this work to track the tensile strength, in practice, due to low-frequency measurement availability from the SOTAX AT4 tablet tester, sensor fusion methods might be required to integrate and efficiently utilize all available plant data. This improved strategy would also require additional studies to determine how frequently to collect measurement data from the SOTAX AT4 tablet tester, due to the destructive nature of the testing method.

**Author Contributions:** Conceptualization, Y.-S.H. and M.Z.S.; methodology, Y.-S.H. and M.Z.S.; investigation: Y.-S.H. and S.B.; writing—original draft preparation, Y.-S.H., M.Z.S. and S.B.; writing— review and editing, M.G., Z.K.N., G.V.R.; funding acquisition, G.V.R., Z.K.N., M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by United States Food and Drug Administration gran<sup>t</sup> number 1U01FD006487-01.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
