Materials

A two-component, direct compression compatible model system was selected containing an API and a tableting excipient. The API was crystalline acetylsalicylic acid (ASA) (particle size: d10 = 133.16 µm; d50 = 532.86 µm; d90 = 1223.77 µm) marketed as Rhodine 3040 (Novacyl, France). The tableting excipient was microcrystalline cellulose (MCC) marketed as Vivapur 200 (particle size: d10 = 113.88 µm; d50 = 247.93 µm; d90 = 449.58 µm), bought from JRS Pharma GmbH (Germany). The two components were mixed in a ratio of 20% ASA to 80% MCC, but the ratio was changed in some measurements to provide disturbance for the system.

Continuous Powder Blending Setup

The investigated powder blending is the same process that we used in the paper of Gyürkés et al. [23] (Fig. 1A). In continuous lines, the ingredients are simultaneously and continuously fed into a blender to ensure proper blend uniformity for the subsequent operations. The MCC was fed with a twin-screw gravimetric feeder (Brabender loss-in-weight feeder, DDW-MD0-MT-1.5 HYD with Congrav OP 1 T operator interface). The screws have high capacity and oval geometry (9.0 mm × 11.5 mm).

Fig. 1

A Powder blending setup and B screw configuration

The ASA was fed with a single-screw feeder (FPS Pharma, Fiorenzuola d’Arda, Italy), which has a round geometry and a lower feeding capacity. It is ideal for low mass flow rates. The feeder has no integrated scale or control system, but the rotation speed can be controlled via a serial connection.

The ASA and MCC were continuously blended in a TS16 QuickExtruder (Quick 2000 Ltd., Hungary) twin-screw multipurpose equipment using a 25 L/D ratio screw with 16 mm diameter. The equipment proved feasible for powder blending [19, 23, 51]. The screw contained conveying and kneading elements (Fig. 1B), but any changes in particle size and flow dynamics were not experienced. The excipients were fed with different feeders through the same feeding chute.

After the twin-screw blender, the product was collected on a conveyor belt, which was operated at 18 cm/min speed. The NIR probe was mounted on the conveyor belt. Generally, the probe was placed as close to the extruder die to achieve real-time results. However, in some measurements, the probe was placed 15 cm from the die modeling process with a longer dead time. The configuration, except the NIR probe in the aforementioned measurements, remained unchanged between the RTD and control measurements.

The prepared blends are simple formulations. More components may require further feeders or multiple coupled powder blending steps for more complicated formulations, which may also make the NIR models complicated.

Feeder Profile

MCC was fed into the blender with the gravimetric loss-in-weight (LIW) feeder, which was controlled with the internal algorithm of the feeder. The feeding rate was measured and shown on a control panel. The feeder has an inbuilt calibration based on a lower and higher feeding rate.

However, the ASA was fed with a volumetric feeder, which has no inbuilt weight measurement but can be controlled externally. We have measured the feeder profile to use as a base for the control. The profile was measured offline, with a gain-in-weight (GIW) setup, and in-line with the continuous powder blending setup, measured with the NIR probe.

The GIW setup consisted of the feeder and an analytical scale (Sartorius L420P). The analytical scale was placed under the feeder, catching the fed ASA. The feeder was controlled externally, setting the speed of the feeder uniformly, reading the weight after the feeder was set at the controlled speed, and reading the GIW every minute for 3 min. The intermediate settings were calculated with the linear fit between two measurements.

NIR Spectrometry and Data Analysis

NIR spectrometry was utilized for the characterization of the process and as the input sensor in the feedback control loop. Bruker MPA FT-NIR spectrometer (Bruker Optik GmbH, Germany) was equipped with a Solvias fiber optic probe (Solvias AG, Switzerland) in reflection mode. NIR spectra were collected with 8 cm−1 resolution at the range of 4000–12,500 cm−1. In the case of calibration, 16 scans were accumulated for each measurement, while in the case of real-time measurements, the spectra were collected without accumulation; therefore, the spectra collection was less than 4 s. The concentration was calculated from the spectra with the PLS calibration used by Gyürkés et al. [23].

The measured concentration curves were smoothed for visualization and quantification purposes but were not used for control. Smoothing was achieved with Savitzky-Golay smoothing based on a 15-data point window and second polynomial order. The controller used the raw data, but the performance of the controller is presented based on the smoothed curve. Smoothing was necessary to reduce the imperfections caused by the fast acquisition time, which derives not only from the inaccuracy of the NIR evaluation but also from the small spot size and, therefore, not representative sample size caused by the NIR probe.

Spectra were collected by OPUS 7.5 software (Bruker). The spectra were evaluated by MATLAB R2020b (9.9, Mathworks, USA) and PLS Toolbox 8.7.1. (Eigenvector Research, USA). The concentration monitoring and control were carried out by in-house developed MATLAB graphical user interfaces (GUI). The GUI runs in real-time, reading the spectra from files, then calculating the concentration with predefined partial least squares (PLS) models. The GUI calculates the control action based on PID control loop and in the case of Smith-predictor based on predefined RTD models and, finally, controls the feeder through a serial connection.

Determination of Residence Time Distribution

The RTD of the continuous powder blender setup was measured with impulse disturbances. MCC was continuously fed to the hopper with gravimetric feeding. Pure ASA was manually added to the system directly into the hooper. To minimize the error generated by the added ASA, the impulse size was set to 0.2 g.

The blend product was measured spectroscopically in real time, and the response function was measured with the GUI by starting a timer at the time of the disturbance and collecting the calculated concentrations.

The response function corresponds to the probability distribution function (PDF) of the RTD after baseline correction and normalization, which was done by dividing the response with the area under curve (AUC) (Eq. 1). The main parameters of the RTD were calculated from the moments (Eqs. 2–4) of the response function calculated with numerical integrations.

$$E\left(t\right)=\frac_^\left(c\left(t\right)-c\left(0\right)\right)^}},$$

(1)

$$\mathrm=_^\mathrm\cdot \mathrm\left(\mathrm\right)\cdot \mathrm,$$

(2)

$$\mathrm=_^-\mathrm\right)}^\cdot \mathrm\left(\mathrm\right)\cdot \mathrm,$$

(3)

where E(t) is the probability distribution function (PDF) of the process in t as time, c(t) is the measured concentration in t, MRT is the mean residence time of the PDF, Var is the variance of the PDF, and σ is the standard deviation of the PDF.

Residence Time Distribution Modeling

The modeling was done with the least squares method after preprocessing of the measured responses.

The preprocessing contained a baseline correction to adjust the bias of the PLS model at zero concentration. The baseline correction is also needed if the concentration before and after the impulse disturbance is not zero, i.e., when the disturbance is done with a pre-blend or during normal operation, with all ingredients fed simultaneously. After the baseline correction, the responses were normalized by numerical integration. The preprocessing was done according to Eq. (1) and directly resulted in the measured PDF.

Tank-in-series (TIS) models were used for RTD modeling. The model is described in Eqs. (5 and 6). Plug flow attribute was described with a dead time (\(_)\). The TIS-like behavior was described with two additional parameters using equal size tanks, the number of the tank (N), and the MRT of the total system (\(_\)).

$$_\left(t\right)=\frac^_}_}}\cdot _}_}\right)}^}_\cdot \left(N-1\right)!},$$

(5)

where \(_\) is the PDF based on the TIS model and \(_\) is the MRT of one individual Tank.

The least squares method was accomplished by finding the minimum of the sum of squared error (SSE), calculated based on Eq. (7), by changing the value of the three parameters of the TIS model.

$$SSE= \sum_^_(t)\right|}__}-E(i)\right]}^,$$

(7)

where i refers to the ith measurement at ti time.

Feedback Control

Two control structures were compared, the simple PI(D) control and the advanced control structure of the Smith predictor. Figure 2 presents the PI(D) control structure, in which the controller calculated the actual feeder speed based on Eqs. (8 and 9).

$$Er\left(t\right)=SP\left(t\right)-c\left(t\right),$$

(8)

$$F\left(t\right)=B\left(SP\left(t\right)\right)+_\cdot Er\left(t\right)+ _\cdot _^Er\left(t\right)dt+_\cdot \frac,$$

(9)

where Er(t) is the error in t time, SP(t) is the setpoint at t time, and c(t) is the concentration at t time calculated from the NIR spectrum. The controller calculates F(t) as the feeder speed at t, from B(SP(t)) as the base value, calculated from the actual setpoint based on the feeder profile, and \(_\), \(_\), and \(_\) are the gains of the P, I, and D components of the PID controller.

Fig. 2

The control structure of the Smith predictor has depicted in Fig. 3. The control is the same PID control as in the previous case, but a complex model computed the error values as follows. The model predicted the concentration at the time of the measurement based on the collected feeder speed. This predicted concentration corresponds to earlier feeder speeds due to the effect of RTD. However, the effect of recent feeder speeds is also predicted with the same model, but with the elimination of the dead time. The input of the controller is calculated based on the model error (measured concentration—model predicted concentration), added to the predicted effect of the actual feeder speeds, and subtracted from the setpoint. The calculation is described in Eqs. (10–13).

$$_\left(t\right)=_\left(F\left(t\right)\right)\times E\left(^}\right),$$

(10)

$$E_\left(t\right) = c\left(t\right) - _\left(t\right),$$

(11)

$$_\left(t\right) = _\left(t - _\right) + E_\left(t\right),$$

(12)

$$Er\left(t\right) = SP\left(t\right)-_\left(t\right),$$

(13)

where \(_\left(t\right)\) is the calculated concentration based on the feeder input (\(F\left(t\right)\)), corrected with the feeder profile (\(_(x)\)) and convolved with the PDF of the measured RTD (\(E\left(^}\right)\)). \(E_\left(t\right)\) is the model error. \(_\left(t\right)\) is the Smith-predicted concentration, which is the model concentration without dead time, corrected with the model error. Finally, the PID control is the same as in Eq. (9), using the error calculated with Eq. (13).

Fig. 3

Smith predictor control structure

Control Tuning

Tuning a control loop is an intensive and important part of implementing feedback control. However, the tuned parameters have to be selected according to the system, the limits of the system, and the anticipated deviations. Some systems may be more sensitive for higher but shorter impulses and tolerate long-term, low amplitude deviation better, whereas the opposite may be true for other systems. Therefore, in some cases, slower but robust control tuning is more feasible than faster, aggressive setups [63]. The Smith predictor has been criticized lately because the control structure is more complicated due to the model, while the improved performance can be marginal [63, 64]. However, it is an applicable alternative if PID fails and provides a clear and understandable control when it is combined with the process dynamics.

Intensive control tuning can be achieved by multiple methods, using labor-intensive trial-and-error method, simple rule based methods such as Ziegler-Nichols tuning. The most accurate methods are optimization-based methods, using mathematical optimization, but these are computationally intensive algorithms and require reliable process models. Nevertheless, the tuned parameters require validation and further adjustments.

In the scope of this publication, we were interested in the feasibility of the RTD model-based Smith predictor and the relation between the two control setups. Therefore, a simple tuning method was used based on the simplified open-loop transfer function. The tuning was feasible as the plant does not contain integrating parts.

CDF of the RTD was used as the open-loop transfer function of the system. The function was simplified to a first-order system, by fitting a linear model on the CDF, which was described by the simplified model dead time (tDSM) and the simplified model time constant (tCSM). The fitting is presented in Fig. 4. The control parameters were calculated from two time constants based on Table 1 [65]. For the Smith predictor, the same tuning was used, except that the dead time of the controller was subtracted from the time constants. This method makes the dead time of the Smith predictor an additional parameter of the control structure, which determines the control gains.

Fig. 4

Fit and parameters of the simplified model, used for control tuning

Table 1 Recommended control gains based on the simple tuning