Model predictive control for continuous pharmaceutical feeding blending units

Journal Pre-proof Model predictive control for continuous pharmaceutical feeding blending units Selma Celikovic, Martin Kirchengast, Jakob Rehrl, Julia Kruisz, Stephan Sacher, Johannes Khinast, Martin Horn

PII:

S0263-8762(19)30561-1

DOI:

https://doi.org/10.1016/j.cherd.2019.11.032

Reference:

CHERD 3917

To appear in:

Chemical Engineering Research and Design

Received Date:

6 September 2019

Revised Date:

8 November 2019

Accepted Date:

25 November 2019

Please cite this article as: Celikovic S, Kirchengast M, Rehrl J, Kruisz J, Sacher S, Khinast J, Horn M, Model predictive control for continuous pharmaceutical feeding blending units, Chemical Engineering Research and Design (2019), doi: https://doi.org/10.1016/j.cherd.2019.11.032

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Model Predictive Control for Continuous Pharmaceutical Feeding Blending Units

Selma Celikovica, Martin Kirchengasta, Jakob Rehrlb, Julia Kruiszb, Stephan Sacherb, Johannes Khinastb,c and Martin Horna

a

Institute of Automation and Control, Graz University of Technology, Graz, Austria

[email protected] [email protected] [email protected] b

Research Center Pharmaceutical Engineering, Graz, Austria

[email protected] [email protected] [email protected] [email protected] c

Institute of Process and Particle Engineering, Graz University of Technology, Graz, Austria

ro of

Corresponding author: Selma Celikovic; e-mail: [email protected]; address: Institute of Automation and Control, Graz University of Technology, Graz, Austria

na

lP

re

-p

Graphical Abstract

  

Several controller strategies for a feeding-blending unit are designed and applied in a simulation and on a real plant Standard controller concepts are compared to a model predictive controller Brief feeder fluctuations are successfully attenuated by adjusting the blender speed Application of a model predictive controller shows advantageous behaviour compared to the other controller strategies

Jo



ur

Highlights

Model Predictive Control for Continuous Pharmaceutical Feeding Blending Units Abstract This paper outlines several control strategies for a feeding blending unit in a continuous pharmaceutical manufacturing setup, with the goal of damping feeder fluctuations by adjusting the blender speed. An advanced model predictive control concept and two standard concepts, PI- and model-based feed-forward control, are described and compared with respect to their performance in the simulation. The concepts which showed satisfactory behaviour in the

simulation are further tested on the real system. The paper also provides a detailed explanation of each implementation step in a real system, i.e., measurement device setup, signal processing, parameter identification and controller design. The properties of control concepts are compared and analysed, particularly in view of further work on the complete direct compaction line.

Keywords: Model predictive control; Continuous pharmaceutical manufacturing; Feeding and blending; Advanced process control; Process modelling

1. Introduction

-p

ro of

Continuous pharmaceutical manufacturing has become a focus of current research in the pharmaceutical industry due to the benefits it brings (Lee et al. 2015). This type of production has many advantages over conventional batch-based production lines, such as higher productivity, lower costs, great flexibility, realtime production and constant product quality. A change from batch to automated continuous production has been promoted in recent years. A continuous plant typically has a centralized control platform, capturing real-time process data and executing control strategies in order to obtain the desired product quality (Singh 2018). In continuous plants, it is crucial that the quality of mixed materials remains the same over time and that the low-quality materials are discharged immediately. This requires control strategies ensuring that the critical process parameters stay within the desired reference values, sustaining the appropriate operation of a continuous production (Rehrl et al. 2016). Providing quality control and risk management, a robust and reliable controller is of great importance for a continuous plant (Su, Moreno, et al. 2018).

lP

re

A continuous manufacturing setup contains diverse devices, such as a feeder, a blender, a tablet-press, etc. An integral component of almost all continuous production lines is a feeding-blending unit (FBU). An FBU consists of two or more feeders and one blender. Its performance has a significant impact on the final product quality. On common continuous tableting lines, a tablet press is placed below the FBU and supplied with a specified, constant inlet mass flow. In this configuration, the inlet mass flow of the tablet press corresponds to the blender outlet mass flow. Since tablet press disturbances have a significant influence on the quality attributes of the final product (Su, Bommireddy, et al. 2018), the inlet mass flow of the tablet press should remain as close to the desired reference value as possible.

Jo

ur

na

In this work, the emphasis was placed on the design and implementation of a control strategy for FBUs, focusing on control for the blender outlet mass flow. The device setup consists of an FBU with an unusual configuration (one feeder and one blender) and a scale at the bottom of the unit. A single feeder only is sufficient since the blender outlet mass flow is the only parameter of interest, while the content uniformity and concentration of the blend are not analysed. Content uniformity and concentration are of great importance, and it would be logical to extend the proposed concept by considering those properties at a later stage. The effects of feeder fluctuations, which cause temporary deviations in the blender inlet mass flow and a deviation of the blender outlet mass flow from the reference value, should be compensated for in order to ensure the desired mass flow at the blender outlet. The effect of these fluctuations can be eliminated by adjusting the blender speed. Several control strategies were designed, implemented and compared with respect to their performance at a real plant. The system behaviour is characterized mainly by the blender dynamics. There are several approaches to derive the blender model and perform its parameter identification. One possibility is a mathematical model based on physical principles describing the hold-up (mass of material) in the blender (Rehrl et al. 2016). Another one is a direct identification of linear models based on the laboratory measurements, i.e., identifying the parameters via numerical minimization of the model error with regard to the measurement data. The resulting model with linear behaviour is valid in the domain around the chosen operating point. In our work, this model is used for simulation and parameterization of three control concepts suggested and simulated in (Rehrl et al. 2016). The first-mentioned control concept is a proportional-integral (PI) controller with an additional feedback loop. A standard PI controller does not consider the effect of the blender speed change on the hold-up in the blender, e.g., a significant deviation of the blender speed

from the operating point could lead to an undesired change in the hold-up and the mixing regime. To that end, an additional feedback loop extends the PI controller (Section 4.1). The second possibility is a model-based feed-forward controller (Section 4.2), which uses the measured feeder disturbance value and the model to compute the required change in the blender speed. The advantage of this is that the outlet mass flow of the blender does not have to be measured, which would be relatively difficult to accomplish on a continuous production line. Similar to the PI controller, however, this concept does not consider the effect of the blender speed on the hold-up. The third approach is a model predictive controller (MPC), (Section 4.3). MPC has many advantages over the other two standard controller concepts, the main one being the ability to keep the process parameters within the set boundaries. Bhaskar, Barros and Singh, 2017 report that its application results in a better closed-loop process performance, despite the MPC requiring a higher computational effort than standard controllers.

Jo

ur

na

lP

re

-p

ro of

The remainder of this paper is structured as follows: Section 2 contains a detailed explanation of the individual components and the device setup. In Section 3, transfer functions and a state-space representation of the blender model are presented. The signal processing steps are explained and the hold-up term is introduced. This Section also contains a comparison between the simulated model and real plant behaviour. In Section 4, a PI controller with a feedback loop, a model-based feed-forward controller and an MPC are introduced. These control concepts are compared to each other in the simulation. The concepts which perform well in the simulation are further tested at a real plant. Finally, Section 5 concludes the paper, providing a summary of our work and outlining possible future research.

Figure 1: Experimental device setup with a feeding blending unit and a scale at the real plant (left), sketch of the device setup (right). 2. Device setup Figure 1 shows the device setup, which consists of a feeder, K-Tron K-CL-KT20 with K-SFS-24 scale (“Coperion GmbH” 2019), a blender, Hosokawa Modulo-mix (“Hosokawa GmbH” 2019) and a custom-built

catch scale, HBM PW22C3 single point load cell (“HBM GmbH” 2019). The feeder operates in the gravimetric mode, supplying the blender with desired inlet mass flow 𝑚̇𝑖𝑛 . The fed material is processed in the blender at a specified rotational speed 𝑛 and subsequently falls on the catch scale. The blender’s outlet mass flow 𝑚̇𝑜𝑢𝑡 is indirectly measured via the catch scale. As shown in Figure 1, the device setup comprises only one feeder and one blender, since the purpose of this study was to examine the blender’s outlet mass flow and hold-up behaviour without analysing the content uniformity and concentration of the final product. Nominal values of the FBU parameters (values in steady state), 𝑚̇𝑖𝑛,𝑛𝑜𝑚 = 5𝑘𝑔/ℎ, 𝑚̇𝑜𝑢𝑡,𝑛𝑜𝑚 = 5𝑘𝑔/ℎ and 𝑛𝑛𝑜𝑚 = 800𝑟𝑝𝑚, are empirically selected via laboratory tests. Pre-blend material used in laboratory measurements contains vitamin D3 (5.5%), mannitol (60.7%), microcrystalline cellulose (4.4%), silica colloidal anhydrous (0.9%), maize starch fine powder (16.2%), sodium starch glycolate (6.4%), talk (5%) and Mg stearate (0.9%). It should be noted that after a material change, model identification and controller parameter design have to be adjusted via new laboratory measurements.

lP

re

-p

ro of

The hardware configuration and the device control were implemented using XAMControl Automation Software (“Evon GmbH” 2019). A real-time database allows updating and saving all data of interest during the test run measurements applying the user-defined sampling time. The saved variables are available at any time and can be used in subsequent data processing. The feeder is connected to XAMControl via a PBpro ETH remote Profibus interface (“Softing GmbH” 2019). The catch scale and the blender are connected via a B&R Bus Controller and analog I/O modules (“B&R Industrial Automation” 2019).

3. Blender model

na

Figure 2: Blender model with 2 inputs and 1 output.

Jo

ur

This Section describes the identification of a continuously-operated blender model. One possible approach is deriving a mathematical model based on the physical principles that describe the mass of material in the blender. The result, however, is a non-linear model (Rehrl et al. 2016), while the purpose of this work is a controller design for specified operating points, which would require linearization of the non-linear model in these points. To address this issue, a direct model identification based on the laboratory measurements made in the experimental device setup was chosen. Parameters are identified through numerical minimization of a model error with regard to the measurement data. Figure 2 shows the block diagram of the blender model with two inputs: blender speed 𝑛 as an actuating signal and inlet mass flow 𝑚̇𝑖𝑛 as a measured disturbance. The symbol Δ is used to indicate deviations from the operating (equilibrium) point, which is denoted by the subscript “eq”. The FBU nominal values introduced in Section 2 are used as equilibrium points for model identification. 3.1 Signal Processing The blender’s outlet mass flow 𝑚̇𝑜𝑢𝑡 is not directly measurable on the catch scale. The measured mass of the material falling on the scale (𝑚𝑜𝑢𝑡 ) and time difference between the samples ∆𝑡 (XAM Control sampling time, ∆𝑡 = 𝑇𝑠 = 0.35s) can be applied to compute the outlet mass flow as:

𝑚̇𝑜𝑢𝑡 (𝑡) =

𝑑𝑚𝑜𝑢𝑡 (𝑡) 𝑚𝑜𝑢𝑡 (𝑡) − 𝑚𝑜𝑢𝑡 (𝑡 − ∆𝑡) ≈ 𝑑𝑡 ∆𝑡

(1)

Due to a high sampling rate and the noise-amplifying nature of the derivative operator, this signal is not suitable for any further work. However, the outlet mass flow signal is crucial for model identification and controller design. As such, an adequate filter must be designed with low-pass filtering attenuating high frequency noise. This type of filter design is termed ‘frequency selective filtering’.

na

lP

re

-p

ro of

The type and parameters of the filter are summarized in Table 1. A digital low-pass finite impulse response (FIR) filter was used for online real-time measurements. Filter order 𝑛𝑓𝑖𝑙𝑡 was chosen as a compromise between rejecting the undesired high-frequency noise and avoiding a phase delay in the filtered signal. The cut-off frequency 𝑓𝑐𝑢𝑡−𝑜𝑓𝑓 in Table 1 was selected experimentally. The PI controller and MPC observer introduced in Section 4 use this filtered signal as a measured process variable. Forward-backward filtering was applied for model identification purposes since it does not shift the signal while filtering and causes no time delay in the resulting models. Yet this type of filtering is non-causal and can be implemented only in offline processing. Figure 3 shows the comparison between unfiltered and filtered (in two ways) outlet mass flows. It is evident that the unfiltered signal (Equation 1) cannot be used. Signals obtained using the second filter are smoother and introduce no time delay, as opposed to the signal obtained by applying the first filter.

ur

Figure 3: Comparison between unfiltered, FIR filtered (filter 1, Table 1) and forward-backward filtered (filter 2, Table 1) outlet mass flows via inlet mass flow change and constant blender speed (left), zoom of the marked process event (right). 3.2 Plant transfer functions

Jo

It is assumed that the plant behaviour is linear in the case of small deviations from the steady-state operating points. The transfer functions 𝑃(𝑠) and 𝑃𝑑 (𝑠) are developed on this basis in this Section. The structure of the transfer functions is defined and the parameters are identified based on the laboratory experiments. Since transfer functions 𝑃(𝑠) and 𝑃𝑑 (𝑠) are identified using the forward-backward filtered outlet mass flow, there is no time delay. 𝑃𝑑𝑒𝑙𝑎𝑦 (𝑠) and 𝑃𝑑,𝑑𝑒𝑙𝑎𝑦 (𝑠) are established directly using the real-time filtered data for filter 1 in Table 1. Time-delay-free transfer functions are used to design the disturbance feed-forward controller concept (described in Section 4.2), which needs no feedback of measured outlet mass flow. 𝑃𝑑𝑒𝑙𝑎𝑦 (𝑠) and 𝑃𝑑,𝑑𝑒𝑙𝑎𝑦 (𝑠) are applied in the PI controller concept (given in Section 4.1), which requires online filtered outlet mass flow as a measured process variable. Transfer function 𝑃(𝑠) describes the effect of blender speed change ∆𝑛 on the outlet mass flow change ∆𝑚̇𝑜𝑢𝑡 . In the first laboratory experiment, the blender’s inlet mass flow was set to mass flow equilibrium input 𝑚̇𝑖𝑛,𝑒𝑞 while the blender speed was varied around the equilibrium input 𝑛𝑒𝑞 . The outlet mass flow 𝑚̇𝑜𝑢𝑡 is determined as described in Section 3.1. At t=0s the blender speed and inlet mass flow are set to nominal

values. At t=285s the blender speed is increased by 200rpm, and at t=765s it is set to the nominal value. At t=1245s the blender speed is then decreased by 200rpm, and at t=1725s it is set back to nominal value. The difference between the measured blender speed 𝑛 and the equilibrium input 𝑛𝑒𝑞 was used as an input and the difference between outlet mass flow 𝑚̇𝑜𝑢𝑡 and equilibrium output 𝑚̇𝑜𝑢𝑡,𝑒𝑞 was an output signal for the System Identification Tool in Matlab (“MathWorks GmbH” 2019). The order of identified transfer function was selected experimentally to be a compromise between the system’s accuracy and complexity. The outcome was a second order bi-proper transfer function. Since the system properties are known, the common mode gain of the transfer function must be zero. To that end, we set 𝑏0 = 0 to ensure that both transfer functions have a zero at s = 0. This means that it is possible to obtain the same outlet mass flow at various blender speeds (in other words, the blender speed has no effect on the outlet mass flow when the FBU operates in a steady state). 𝑃(𝑠) =

̅̅̅̅̅̅̅̅ ∆𝑚 ̇ 𝑜𝑢𝑡 (𝑠) 𝑏2 𝑠 2 + 𝑏1 𝑠 + 𝑏0 𝑏2 𝑠 2 + 𝑏1 𝑠 = 2 = 2 ̅̅̅̅ 𝑠 + 𝑎1 𝑠 + 𝑎0 𝑠 + 𝑎1 𝑠 + 𝑎0 ∆𝑛(𝑠)

𝑃𝑑𝑒𝑙𝑎𝑦 (𝑠) =

(2)

̅̅̅̅̅̅̅̅ ∆𝑚 ̇ 𝑜𝑢𝑡 (𝑠) 𝑏2 𝑠 2 + 𝑏1 𝑠 + 𝑏0 𝑏2 𝑠 2 + 𝑏1 𝑠 = 𝑒 −𝑇𝑑 𝑠 2 = 𝑒 −𝑇𝑑 𝑠 2 ̅̅̅̅(𝑠) 𝑠 + 𝑎1 𝑠 + 𝑎0 𝑠 + 𝑎1 𝑠 + 𝑎0 ∆𝑛

(3)

̅̅̅̅̅̅̅̅ ∆𝑚 ̇ 𝑜𝑢𝑡 (𝑠) 𝑏0 = 2 ̅̅̅̅̅̅̅ 𝑠 + 𝑎1 𝑠 + 𝑎0 ∆𝑚̇𝑖𝑛 (𝑠)

̅̅̅̅̅̅̅̅ ∆𝑚 ̇ 𝑜𝑢𝑡 (𝑠) 𝑏0 = 𝑒 −𝑇𝑑 𝑠 2 ̅̅̅̅̅̅̅ 𝑠 + 𝑎1 𝑠 + 𝑎0 ∆𝑚̇𝑖𝑛 (𝑠)

na

𝑃𝑑,𝑑𝑒𝑙𝑎𝑦 (𝑠) =

lP

𝑃𝑑 (𝑠) =

re

-p

ro of

Transfer function 𝑃𝑑 (𝑠) describes the effect of inlet mass flow change ∆𝑚̇𝑖𝑛 on outlet mass flow ∆𝑚̇𝑜𝑢𝑡 . In the second laboratory experiment, the blender speed was set to blender speed equilibrium input 𝑛𝑒𝑞 and the inlet mass flow was varied around inlet mass flow equilibrium input 𝑚̇𝑖𝑛,𝑒𝑞 . The outlet mass flow 𝑚̇𝑜𝑢𝑡 was computed and filtered using the methods described in Section 3.1. At t=0s the blender speed and the inlet mass flow were set to their nominal values. At t=285s the inlet mass flow was increased by 1.5kg/h, and at t=620s it was decreased by 3kg/h. At t=925s it was increased by 4.5kg/h, and at t=1225s it was set back to the nominal value. The difference between measured inlet mass flow 𝑚̇𝑖𝑛 and equilibrium input 𝑚̇𝑖𝑛,𝑒𝑞 was used as an input, while the difference between measured outlet mass flow 𝑚̇𝑜𝑢𝑡 and equilibrium output 𝑚̇𝑜𝑢𝑡,𝑒𝑞 was an output signal for the System Identification Tool in Matlab. The order of 𝑃𝑑 (𝑠) was determined in the same manner as the transfer function 𝑃(𝑠), i.e., as a compromise between accuracy and complexity. The outcome was a second-order transfer function with low-pass characteristics containing no zeros and two poles. (4) (5)

Jo

ur

Figures 4 and 5 indicate very good conformity of the system responses for estimated transfer functions 𝑃(𝑠) and 𝑃𝑑 (𝑠) with the filtered measurement data (filter 2 from Table 1) from the test runs.

na

lP

re

-p

ro of

Figure 4: Test run 1, change in the blender speed via a constant inlet mass flow and a comparison between the measured and estimated outlet mass flows (left), zoom of the marked process event (right).

3.3 Hold-up

ur

Figure 5: Test run 2, change of inlet mass flow via a constant blender speed and a comparison between the measured and estimated outlet mass flows (left), zoom of the marked process event (right).

Jo

During the test runs the material fed by the feeder can either stay in the blender or fall on the catch scale. The mass of the material inside the blender from this point on referred to as hold-up 𝑚ℎ𝑢 is the difference between the mass fed from the feeder (𝑚𝑖𝑛 ) and the total mass measured via catch scale (𝑚𝑜𝑢𝑡 ) and can be computed as follows: ∆𝑚𝑖𝑛 = 𝑚𝑖𝑛 (𝑡𝑐𝑢𝑟𝑟𝑒𝑛𝑡 ) − 𝑚𝑖𝑛 (𝑡 = 0)

(6.1)

∆𝑚𝑜𝑢𝑡 = 𝑚𝑜𝑢𝑡 (𝑡𝑐𝑢𝑟𝑟𝑒𝑛𝑡 ) − 𝑚𝑜𝑢𝑡 (𝑡 = 0)

(6.2)

𝑚ℎ𝑢 = ∆𝑚𝑖𝑛 − ∆𝑚𝑜𝑢𝑡

(6.3)

Figure 6 shows the hold-up signal during test runs 1 and 2. The signal evolution indicates that, although the plant has been operating in a steady state for a long period of time, the hold-up constantly increases. Due to mass conservation, the hold-up should stay constant if the inlet and outlet mass flow rates are the same. This deviation from the expected behaviour is a consequence of the dust extraction system which is brought into the device setup because of the high level of dust formation (Figure 1). The mass of extracted material and the total mass measured via the catch scale should be added up and considered in the hold-up computation.

Since the extracted mass is not measurable, we make the assumption that its amount has a linear increase over time. Consequently, a slope 𝑘ℎ𝑢 of a linear function shown in Figure 6 is computed and the hold-up value is corrected as follows: 𝑚ℎ𝑢,𝑐𝑜𝑟𝑟 = 𝑚ℎ𝑢 − 𝑘ℎ𝑢 ∙ 𝑡

(7)

Equilibrium value 𝑚ℎ𝑢,𝑒𝑞 is a mean value of corrected value 𝑚ℎ𝑢,𝑐𝑜𝑟𝑟 during the intervals at nominal blender speed. Note: From now on the corrected value of hold-up is referred to as hold-up 𝑚ℎ𝑢 . 𝑚ℎ𝑢,𝑒𝑞 = 58𝑔

(8)

lP

re

-p

ro of

Hold-up depends strongly on the blender speed. A decrease in the blender speed causes an immediate increase in the hold-up and vice versa. A permanent increase or decrease in the blender speed may result in a significant hold-up deviation from the operating point, changing the mixing regime.

Jo

ur

na

Figure 6: Change in the hold-up during test runs 1 (left) and 2 (right).

Figure 7: Blender model with 2 inputs and 2 outputs. 3.4 State space representation This Section introduces a linear state space representation of the FBU, which is used to design a linear MPC. It is based on the measurement data from two laboratory experiments described in the previous Section. Figure 7 shows a sketch of the blender model with two input and two output signals.

Note: Up to this point the outlet mass flow has been considered as the only output of blender model. As mentioned in the previous Section, a change in the blender speed causes an immediate change of the holdup. This system property has a great impact on the effectiveness of the controller actions. To that end, the hold-up was used as the second output of the blender model in this Section. The discrete state space model is given by: 𝒙𝑘+1 = 𝑨𝒙𝑘 + 𝑩∆𝑢𝑘 + 𝑩𝒅 ∆𝑑𝑘

(9.1)

∆𝒚𝑘

(9.2)

= 𝑪𝒙𝑘 + 𝑫∆𝑢𝑘 + 𝑫𝒅 ∆𝑑𝑘

𝒙𝑘 ∈ 𝑅 𝑛𝑠𝑦𝑠𝑡 𝑥1 , ∆𝑢𝑘 ∈ 𝑅 𝑚𝑥1 , ∆𝑑𝑘 ∈ 𝑅 𝑙𝑥1 , ∆𝒚𝒌 ∈ 𝑅 𝑝𝑥1

(9.3)

𝑨 ∈ 𝑅 𝑛𝑠𝑦𝑠𝑡 𝑥 𝑛𝑠𝑦𝑠𝑡 , 𝑩 ∈ 𝑅 𝑛𝑠𝑦𝑠𝑡 𝑥 𝑚 , 𝑩𝒅 ∈ 𝑅 𝑛𝑠𝑦𝑠𝑡 𝑥 𝑙 , 𝑪 ∈ 𝑅 𝑝 𝑥 𝑛𝑠𝑦𝑠𝑡 , 𝑫 ∈ 𝑅 𝑝 𝑥 𝑚 , 𝑫𝒅 ∈ 𝑅 𝑝 𝑥 𝑙

(9.4)

Jo

ur

na

lP

re

-p

ro of

Input signal ∆𝑢𝑘 is the difference between measured blender speed 𝑛 and equilibrium input 𝑛𝑒𝑞 (𝑚 = 1, i.e., one actuating signal). The measured disturbance signal ∆𝑑𝑘 is the difference between measured inlet mass flow 𝑚̇𝑖𝑛 and equilibrium input 𝑚̇𝑖𝑛,𝑒𝑞 (𝑙 = 1, i.e., one disturbance signal). The output vector consists of the difference between offline measured outlet mass flow 𝑚̇𝑜𝑢𝑡 (Table 1, filter 2) and equilibrium output 𝑚̇𝑜𝑢𝑡,𝑒𝑞 , and the difference between computed hold-up 𝑚ℎ𝑢 and its equilibrium value 𝑚ℎ𝑢,𝑒𝑞 (𝑝 = 2, i.e., two outputs). State vector 𝒙 contains 𝑛𝑠𝑦𝑠𝑡 states with no physical meaning. The model parameters (i.e., matrices 𝑨, 𝑩, 𝑩𝒅 , 𝑪, 𝑫, 𝑫𝒅 ) are computed in Matlab. The estimation is based on minimizing the error between the measurement and the model for the given signals. The order is experimentally selected during the estimation, and a state space model of third order (i.e., 𝑛𝑠𝑦𝑠𝑡 = 3) is chosen as a compromise between the system’s accuracy and complexity. The state space model parameters are given in Table 3. Figures 8 and 9 show a very good agreement between the system’s response for the estimated state space model and the filtered measurement data from test runs. Note: In order to keep the notation for the MPC design simple, the delta symbols are omitted hereinafter.

Figure 8: Test run 1, change in the blender speed at a constant inlet mass flow and a comparison between the measured and state-space-model-estimated outlet mass flows and hold-ups (left), zoom of the marked process event (right).

ro of

Figure 9: Test run 2, change in the inlet’s mass flow at a constant blender speed and a comparison between the measured and state-space-model-estimated outlet mass flows and hold-ups (left), zoom of the marked process event (right).

-p

3.5 Validation test runs

Jo

ur

na

lP

re

In order to verify the efficiency of the transfer functions and state space model identified in previous Sections, two validation test runs were performed on the real plant. In the first validation test run inlet mass flow stays constant during the variation of the blender speed. Measured value of the outlet mass flow is compared to state-space-model and transfer-function-estimated outlet mass flow. Hold-up is computed from measurement data and compared to the state-space-model-estimated hold-up. In the second validation test the same signals of interest, i.e., measured- and estimated outlet mass flow and hold-up, are compared during the inlet mass flow variation at the constant blender speed. In order to check a model efficiency in a wider range both tests are performed close, but not exactly at the nominal blender operating point. Figure 10 indicates very good conformity of the estimated system outputs with the measurement data from the test runs.

Figure 10: Comparison between measured and estimated hold-up and outlet mass flow values during the validation test run 1 (left) and validation test run 2 (right).

4. Controller design The feeder fluctuations cause a deviation between the real inlet mass flow and the desired set point at the blender inlet. If the blender operates at a constant speed, its inlet fluctuations are transferred to its outlet. In order to achieve a constant mass flow at the blender outlet, the effect of these fluctuations should be eliminated. To that end, the control law uses the blender speed as a manipulated variable. This Section presents three control strategies that were examined in the simulation. The strategies which showed satisfactory behaviour in the simulation are further tested at the real plant. Functionality of the controller action is validated via following test scenario. At time t=0s the inlet mass flow is set to the nominal value. At time t=100s an artificially-introduced feeder disturbance occurs and the mass flow increases by 0.5kg/h. At t=150s it decreases to 4.5kg/h, and at t=200s it is reset to the reference value. At this point it should be noted, that the selected disturbances are chosen much higher than typically expected deviations in such a setup. However, if the proposed concepts are able to supress even those types of disturbances, then smaller disturbances which typically occur in real situations will be attenuated as well. Furthermore, these types of mass flow fluctuations will of course also affect the content uniformity of the final blend. However, the current work focuses on mass flow control and the aspect of content uniformity is not investigated. The following four scenarios were evaluated in this test case:

ro of

 

no control action control strategy 1: proportional-integral (PI) controller with feedback loop (i.e., first order lag compensator) control strategy 2: model-based disturbance feed forward control strategy 3: model-predictive control (MPC) with observer

 

-p

4.1 Control strategy 1: PI controller with feedback loop (first order lag compensator)

1 + 𝑇𝑁 𝑠 𝑘𝑢 (1 + 𝑇𝑁 𝑠(1 +

1 )) 𝐾𝑝 𝑘𝑢

ur

𝑅(𝑠) =

na

lP

re

A PI controller with feedback loop is introduced in this Section. As previously mentioned in Section 3.2, due to the structure of the transfer function 𝑃𝑑𝑒𝑙𝑎𝑦 (𝑠), the blender speed has no influence on the outlet mass flow when the FBU operates in a steady state. Nevertheless, the blender speed directly affects the hold-up in the blender (Section 3.3). The PI controller uses the online- filtered outlet mass flow as the measured process variable. Inaccuracies in the outlet mass flow measurements cause a permanent discrepancy between the measured mass flow and the desired set point. The PI controller reacts to this deviation up to the saturation limits. This behaviour leads to a stationary blender speed at one of the saturation limits. Consequently, the hold-up in the blender at steady state is also far away from the desired operating point. To counteract this, an additional feedback loop 𝑘𝑢 is introduced into the controller structure. The purpose of this feedback loop is to keep the blender speed close to the nominal value in a steady state. The resulting controller transfer function 𝑅(𝑠) is given as: (10)

Jo

Unlike the standard PI controller, the new transfer function has no pole in 𝑠 = 0. This means that the controller error cannot be reduced to zero. Choosing parameter 𝑘𝑢 is a compromise between the controller’s error and the steady-state value of the blender speed. Higher values of 𝑘𝑢 cause a lower change in the blender speed and a larger deviation from the outlet mass flow set point. Yet if parameter 𝑘𝑢 is too small a great deviation from the blender speed reference value in a steady state occurs. Figure 11 shows a block diagram of control strategy suggested in this Section. The controller’s parameters are empirically selected via simulation tests (𝐾𝑝 = 7, 𝑇𝑁 = 10, 𝑘𝑢 = 0.01) and the blender speed is kept between 𝑛𝑚𝑖𝑛 = 600 and 𝑛𝑚𝑎𝑥 = 1000 . Figure 12 illustrates the outlet mass flow signal during the artificially-introduced feeder disturbance in the simulation and at the real plant when no control action takes place (left) and the simulated controller action using the control strategy 1 (right). Although the deviation in the outlet mass flow is reduced, the simulation indicates that this controller is not suitable for the measurement configuration in question: it cannot keep the outlet mass flow close enough to the reference value. The reason for this inefficiency is the use of the blender outlet mass flow as a measured process variable. Due to a significant filter time delay in this measurement, the controller-induced change to the blender speed occurs too late and does not yield the desired results. Yet this concept could be used for other materials that are better-suited for the mass flow

measurement (lower dust formation). Due to the unsatisfactory simulation results, this concept was not implemented on a real system.

lP

re

-p

ro of

Figure 11: Block diagram of control strategy 1.

na

Figure 12: No control action in the simulation and at the real plant (left), assessing the controller action using the control strategy 1 in the simulation (right). 4.2 Control strategy 2: disturbance feed-forward

ur

In this Section, a control strategy based on a feed-forward controller is presented. The measurable feeder disturbance ∆𝑚̇𝑖𝑛 is defined as:

Jo

∆𝑚̇𝑖𝑛 = 𝑚̇𝑖𝑛 − 𝑚̇𝑖𝑛,𝑒𝑞

(11.1)

and used to compute the required blender speed variation. The expected change in the outlet mass flow ∆𝑚̇𝑜𝑢𝑡,𝑝𝑟𝑒𝑑 caused by this disturbance is predicted via transfer function 𝑃𝑑 (𝑠) with known input ∆𝑚̇𝑖𝑛 . ∆𝑚̇𝑖𝑛 → 𝑃𝑑 (𝑠) → ∆𝑚̇𝑜𝑢𝑡,𝑝𝑟𝑒𝑑

(11.2)

∆𝑚̇𝑜𝑢𝑡,𝑝𝑟𝑒𝑑 is the deviation from the equilibrium value at the blender outlet. The change in the blender speed necessary for achieving the same deviation at the blender outlet can be computed via inverse transfer function 𝑃(𝑠)−1 . ∆𝑚̇𝑜𝑢𝑡,𝑝𝑟𝑒𝑑 → 𝑃(𝑠)−1 → ∆𝑛

(11.3)

In order to compensate for this deviation, the controller applies the computed blender speed of the opposite sign.

This way, the disturbance feed-forward transfer function is given by: 𝐹(𝑠) = −

̅̅̅̅ ∆𝑛(𝑠) 𝑃𝑑 (𝑠) 59.803(𝑠 2 + 0.335𝑠 + 0.03261) =− = − ̅̅̅̅̅̅̅ (𝑠 + 0.5795)(𝑠 + 0.3765)(𝑠 + 1.309)𝑠 𝑃(𝑠) ∆𝑚̇𝑖𝑛 (𝑠)

(12)

Figure 13 shows a block diagram of the control strategy. The design of 𝐹(𝑠) is based on experimentallyestimated transfer functions. Transfer function 𝑃(𝑠) has a zero in 𝑠 = 0, which results in a pole at 𝑠 = 0 in 𝐹(𝑠). Consequently, the feed-forward transfer function is non BIBO-stable (BIBO: bounded-input, boundedoutput, (Chen 1998) ).To achieve BIBO stability, the pole in 𝑠 = 0 is replaced by a pole of 𝑠 = −0.002. The location of this pole is empirically selected based on the properties of the original transfer function. The new BIBO stable feed-forward transfer function 𝐹𝑏𝑖𝑏𝑜 (𝑠) is computed as: 𝐹𝑏𝑖𝑏𝑜 (𝑠) = −

59.803(𝑠 2 + 0.335𝑠 + 0.03261) (𝑠 + 0.5795)(𝑠 + 0.3765)(𝑠 + 1.309)(𝑠 + 0.002)

(13)

na

lP

re

-p

ro of

Both transfer functions are tested in simulation and at a real plant. 𝐹(𝑠) and 𝐹𝑏𝑖𝑏𝑜 (𝑠) are discretized via Tustin’s method and applied to a real system in a form of difference equations in the discrete time domain. Figure 14 shows the simulated and the real controller action during the artificially-introduced feeder disturbance. The measurement data obtained for the simulated and measured process are in very good agreement. Both controllers (even a BIBO-stable controller with a changed pole) keep the outlet mass flow close to the given reference value during the brief disturbance. Yet the blender speed is quickly driven to the minimal suggested value without considering the blender hold-up.

Jo

ur

Figure 13: Block diagram of control strategy 2.

Figure 14: Assessing the controller’s action using the non-BIBO stable transfer function 𝐹(𝑠) (left) and the BIBO stable transfer function 𝐹𝑏𝑖𝑏𝑜 (𝑠) (right) in the simulation and at the real plant.

4.3 Control strategy 3: model predictive control The disturbance feed-forward and PI controllers described in the previous Sections keep the blender outlet mass flow close to its reference value without considering the blender speed and its effect on the hold-up during the process. As mentioned above, this drawback may have undesired consequences, e.g., a change in the mixing regime. An MPC can address these disadvantages. Of course, the introduced constraints on the blender speed will limit the achievable outlet mass flow control performance. An MPC operates like a good chess player: it plans a number of future moves based on the current conditions and possible future developments, using a dynamic process model to predict the future and optimizing a customized objective function. The entire process can be summarized as follows: The first step is to obtain the current state 𝒙𝑘 either via measurement or from an observer. Next, the controller solves an optimization problem over prediction horizon 𝑛𝑝 starting with the current time step via a numerical minimization algorithm. The control sequence over control horizon 𝑛𝑐 is computed by this means. Subsequently, the first value of this sequence is applied as an actuating signal. Finally, a new system state is obtained (once again via measurement or from an observer) and the optimization is repeated for the next step.

ro of

To design a linear MPC, the following system properties have to be known: 1. A discrete time model of the plant: 𝒙𝑘+1 = 𝑨𝒙𝑘 + 𝑩𝑢𝑘 + 𝑩𝒅 𝑑𝑘 𝒚𝑘

(14.1)

= 𝑪𝒙𝑘 + 𝑫𝑢𝑘 + 𝑫𝒅 𝑑𝑘

(14.2)

re

2. The actuating signal from the last iteration: 𝑢𝑘−1

-p

(Note: Equations 14.1 and 14.2 correspond to Equations 9.1 and 9.2 in Section 3.4, except that ‘∆’ symbols are omitted. The parameters are provided in Table 3.)

3. Current state vector 𝒙𝑘

lP

As already mentioned, the state vector 𝒙𝑘 of the blender model cannot be measured. Therefore, a linear state observer is derived, as shown in the next Section. 4. Reference value 𝒓𝒌 of output signal 𝒚𝒌

na

The reference value of the output signal, i.e., of the outlet mass flow is set to zero in order to assure a constant value at the blender’s outlet.

𝑛𝑝

ur

A control sequence over a control horizon 𝑛𝑐 (𝑢̂𝑘 , 𝑢̂𝑘+1 , … , 𝑢̂𝑘+𝑛𝑐 −1 ), which minimizes a suitable objective function 𝐽, is computed. The formulation of this function is an important step in designing an MPC. Its form, e.g., linear or quadratic, and its parameters directly influence the system’s performance. The chosen quadratic objective function 𝐽 is given as follows: 𝑛𝑐

Jo

𝑇 𝑇 ̂𝑘+𝑖 − 𝒓𝑘+𝑖 )𝑇 𝑸𝒊 (𝒚 ̂𝑘+𝑖 − 𝒓𝑘+𝑖 )] + ∑[𝑢̂𝑘+𝑖−1 𝐽 = ∑[(𝒚 𝑅𝑢,𝑖 𝑢̂𝑘+𝑖−1 + ∆𝑢̂𝑘+𝑖−1 𝑅∆,𝑖 ∆𝑢̂𝑘+𝑖−1 ] 𝑖=1

with:

𝑖=1

̂𝑘+𝑖 𝒚

output at time instant 𝑘 + 𝑖

𝒓𝑘+𝑖

output reference value at time instant 𝑘 + 𝑖

𝑢̂𝑘+𝑖−1

input at time instant 𝑘 + 𝑖 − 1

∆𝑢̂𝑘+𝑖−1

change in input at time instant 𝑘 + 𝑖 − 1, i.e., ∆𝑢̂𝑘+𝑖−1 = 𝑢̂𝑘+𝑖−1 − 𝑢̂𝑘+𝑖−2

𝑸𝒊

output error weighting matrix

𝑅𝑢,𝑖

input (energy) weighting parameter

(15)

𝑅∆,𝑖

input rate weighting parameter

Note: ^ symbol indicates the predicted values by the MPC design. The objective function consists of three terms. The first term relates to the difference between the predicted output and the output reference value, i.e., the control error over the prediction horizon. The second and third terms penalize the actuating signal and a change in it over the control horizon (the same prediction and control horizons are chosen). A choice of 𝑸𝒊 , 𝑅𝑢,𝑖 and 𝑅∆,𝑖 has an impact on weighting the control error, the actuating signal and a change in the actuating signal. By this means a change in the blender speed is introduced in an objective function and tuned via parameter 𝑅∆,𝑖 . A suitable choice of 𝑅∆,𝑖 prevents any undesired or eruptive changes in the blender speed. Hold-up 𝑚ℎ𝑢 is now considered as a second output of the system and its tracking behaviour is controlled via a choice of matrix 𝑸𝒊 . The optimization function given in Equation 15 can be transformed into: ̅ 𝒌 𝑇 𝑯 ∆𝒖 ̅ 𝒌 + 2 𝒇𝑻 ∆𝒖 ̅𝒌 𝐽 = ∆𝒖

(16)

ro of

Matrix 𝑯 and vector 𝒇𝑻 are derived using a recursive representation of system equations and calculated based on the system’s properties {𝑨, 𝑩, 𝑩𝒅 , 𝑪, 𝑫, 𝑫𝒅 , 𝑢𝑘−1 , 𝑥𝑘 , 𝑟𝑘 } and the MPC parameters {𝑛𝑝 , 𝑛𝑐 , 𝑸𝒊 , 𝑅𝑢,𝑖 , 𝑅∆,𝑖 }. The MPC tuning, i.e., choosing the MPC parameters, is performed experimentally by the controller tests at the real plant. The best results for the parameters are shown in Table 4. ̅𝒌: The result of the optimization is a control sequence ∆𝒖

(17)

-p

∆𝑢̂𝑘 ∆𝑢̂𝑘+1 ̅𝒌 = [ ∆𝒖 ] ⋮ ∆𝑢̂𝑘+𝑛𝑐 −1

re

In order to provide the actuating signal 𝑢̂𝑘 to be sent to the plant, the following equation is implemented: 𝑢̂𝑘 = 𝑢𝑘−1 + ∆𝑢̂𝑘

(18)

lP

The greatest advantage of this concept is a relatively simple consideration of constraints. The constraints to the blender speed value and blender’s speed change are given in a form of linear inequalities as: 𝑢𝑚𝑖𝑛 ≤ 𝑢̂𝑘+𝑖 ≤ 𝑢𝑚𝑎𝑥

∀ 𝑖 = {0, … , 𝑛𝑐 − 1}

(19)

∆𝑢𝑚𝑖𝑛 ≤ ∆𝑢̂𝑘+𝑖 ≤ ∆𝑢𝑚𝑎𝑥

∀ 𝑖 = {0, … , 𝑛𝑐 − 1}

(20)

na

̅𝒌: These can be written as the linear constraints of the vector ∆𝒖 ̅𝒌 ≤ 𝒘 𝑾 ∆𝒖

(21)

Jo

ur

Linear constraints can easily be taken into account when minimizing the objective function 𝐽 using state-ofthe-art numerical solvers (Ferreau et al. 2014) to avoid any undesired change in the blender speed that could cause a change in the mixing regime. The desired constraints of blender speed and blender speed change are provided in Table 4. The block diagram of this control strategy with the MPC and the observer is shown in Figure 15.

Figure 15: Block diagram of control strategy 3. 4.3.1 Linear state observer

ro of

One requirement for the MPC design is to know the current state vector 𝒙𝑘 . Since the states of the blender model cannot be measured, an adequate observer has to be derived. The Luenberger state observer provides an estimate of the internal state using the measured values of input and output signals of the real system. A detailed introduction to Luenberger observers can be found in (G.F. Franklin, J.D. Powell 2002). The Luenberger observer for the blender model (Equation 14) reads as: ̂𝑘+1 = 𝑨 𝒙 ̂𝑘 + 𝑩𝑢𝑘 + 𝑩𝒅 𝑑𝑘 + 𝑳(𝑦 𝒙 ⏟ 1𝑚𝑒𝑎𝑠,𝑘 − 𝑦̂1,𝑘 )

-p

𝒄𝒐𝒓𝒓𝒆𝒄𝒕𝒊𝒐𝒏 𝒕𝒆𝒓𝒎

̂𝑘 = 𝑪𝒙 ̂𝑘 + 𝑫𝑢𝑘 + 𝑫𝒅 𝑑𝑘 𝒚

(22.1) (22.2)

re

with an estimation state vector containing 𝑛𝑠𝑦𝑠𝑡 estimated states: 𝑇

̂𝑘 = [𝑥̂1,𝑘 … 𝑥̂ 𝑛𝑠𝑦𝑠𝑡,𝑘 ] 𝒙

(23)

𝑇

̂̇ 𝑜𝑢𝑡,𝑘 ∆𝑚 ̂𝑘 = [𝑦̂1,𝑘 𝑦̂2,𝑘 ] = [ ∆𝑚 𝒚 ̂ ℎ𝑢,𝑘 ]

𝑇

lP

and an estimation output vector consisting of estimated outlet mass flow and estimated hold-up: (24)

̂𝑘 𝒆𝑘 = 𝒙𝑘 − 𝒙

na

Furthermore, the estimation error 𝒆𝑘 is defined as the difference between the actual and estimated state vectors (25)

The error dynamic can be written as:

ur

𝒆𝑘+1 = (𝑨 − 𝑳𝑪𝟏 )𝒆𝑘

(26)

Jo

The observer contains a copy of the system model and a correction term. Measurable system output 𝑦1𝑚𝑒𝑎𝑠,𝑘 = ∆𝑚̇𝑜𝑢𝑡,𝑘 is compared to observer output 𝑦̂1,𝑘 and weighted using design parameter L (second output of system 𝑦2,𝑘 = ∆𝑚ℎ𝑢,𝑘 is not considered in this term since it cannot be measured). By choosing parameter L appropriately, the estimation error 𝒆𝑘 converges to zero. This process is referred to as the eigenvalue placement, i.e., L is chosen such that the poles of error dynamic matrix (𝑨 − 𝑳𝑪𝟏 ) are placed inside the unit circle. The model parameters in (14.1) and (14.2) are determined using the measurement data processed offline (Table 1, filter 2). Therefore, observer output 𝑦̂1,𝑘 considers no delay time or filter impact and yet is compared to real-time measured (and filtered) output 𝑦1𝑚𝑒𝑎𝑠,𝑘 (Table 1, filter 1) with the difference used as a correction term in the state estimation. In order to eliminate this inconsistency, the existing blender model is extended by the filter, i.e., 𝑛𝑓𝑖𝑙𝑡 new states are introduced as described below. The FIR filter in the z-domain reads as: 𝑌𝑓 (𝑧) 𝑌1 (𝑧)

= 𝑓0 + 𝑓1 𝑧 −1 + 𝑓2 𝑧 −2 + ⋯ + 𝑓𝑛𝑓𝑖𝑙𝑡 𝑧 −𝑛𝑓𝑖𝑙𝑡

(27)

where 𝑌1 (𝑧) and 𝑌𝑓 (𝑧) represent z-transformed input and output signals of a filter and 𝑓0 , 𝑓1 , ⋯ , 𝑓𝑛𝑓𝑖𝑙𝑡 the filter coefficients of filter 1 provided in Table 1. In the time domain the following holds: 𝑦𝑓,𝑘 = 𝑓0 𝑦1,𝑘 + 𝑓1 𝑦1,𝑘−1 + 𝑓2 𝑦1,𝑘−2 + ⋯ + 𝑓𝑛𝑓𝑖𝑙𝑡 𝑦1,𝑘−𝑛𝑓𝑖𝑙𝑡

𝑦𝑓,𝑘 = [𝑓0 𝑪𝟏

𝑓1

⋯

𝑓2

𝑓𝑛𝑓𝑖𝑙𝑡 ]

𝒙𝒌 𝑦1𝑑,𝑘 𝑦2𝑑,𝑘 ⋮

(28)

+ 𝑓0 𝐷1 𝑢𝑘 + 𝑓0 𝐷𝑑1 𝑑𝑘

(29)

[𝑦𝑛𝑓𝑖𝑙𝑡𝑑,𝑘 ] In this new notation, the state vector elements [ 𝑦1𝑑,𝑘 , 𝑦2𝑑,𝑘 , ⋯ , 𝑦𝑛𝑓𝑖𝑙𝑡𝑑,𝑘 ] correspond to [𝑦1,𝑘−1 , 𝑦1,𝑘−2 , ⋯ , 𝑦1,𝑘−𝑛𝑓𝑖𝑙𝑡 ] and the parameters 𝑪𝟏 , 𝐷1 and 𝐷𝑑1 to the first row of the parameters 𝑪, 𝑫 and 𝑫𝒅 . The following extended blender model is obtained:

𝑓1

𝑓2

⋯

𝑓𝑛𝑓𝑖𝑙𝑡 ]

𝑪𝒇

𝒙𝒌 𝑦1𝑑,𝑘 𝑦2𝑑,𝑘 ⋮ [𝑦𝑛𝑓𝑖𝑙𝑡𝑑,𝑘 ]

𝒇,𝒅

+ 𝑓⏟ ⏟0 𝐷𝑑1 𝑑𝑘 0 𝐷1 𝑢𝑘 + 𝑓 𝐷𝑓

In a compressed form:

𝐷𝑓,𝑑

lP

𝒙𝑓,𝑘+1 = 𝑨𝒇 𝒙𝑓,𝑘 + 𝑩𝒇 𝑢𝑘 + 𝑩𝒇,𝒅 𝑑𝑘 𝑦𝑓,𝑘 = 𝑪𝒇 𝒙𝑓,𝑘 + 𝐷𝑓 𝑢𝑘 + 𝐷𝑓,𝑑 𝑑𝑘

𝒇

re

𝑦𝑓,𝑘 = ⏟ [𝑓0 𝑪𝟏

𝒙𝑓,𝑘

𝑨𝒇

ro of

𝒙𝑓,𝑘+1

0 0 ⋯ ⋯ 0 𝒙𝑘 𝑩𝒅 𝑩 0 0 ⋯ ⋯ 0 𝑦1𝑑,𝑘 𝐷1 𝐷𝑑1 1 0 ⋯ ⋯ 0 𝑦2𝑑,𝑘 + 0 𝑢𝑘 + 0 𝑑𝑘 0 1 0 ⋯ 0 ⋮ ⋮ ⋮ 0 0 ⋱ ⋱ ⋮ 𝑦 [⏟ ] [ 0 ] ⏟ [ ] 0 𝑛𝑓𝑖𝑙𝑡𝑑,𝑘 ⏟ 0 ⋯ 0 1 0] 𝑩 𝑩

-p

𝑨 𝑪𝟏 = 0 ⋮ ⋮ [𝑦𝑛𝑓𝑖𝑙𝑡𝑑,𝑘+1 ] [ ⏟ ⏟0 𝒙𝑘+1 𝑦1𝑑,𝑘+1 𝑦2𝑑,𝑘+1 ⋮

(30.1)

(30.2)

(31.1) (31.2)

Analogous to (22.1) and (22.2), the Luenberger observer is written as: (32.1)

̂𝑓,𝑘 + 𝐷𝑓 𝑢𝑘 + 𝐷𝑓,𝑑 𝑑𝑘 𝑦̂𝑓,𝑘 = 𝑪𝒇 𝒙

(32.2)

𝑇

ur

̂𝑘 = [𝑦̂𝑓,𝑘 𝑦̂2,𝑘 ] 𝒚

na

̂𝑓,𝑘+1 = 𝑨𝒇 𝒙 ̂𝒇,𝒌 + 𝑩𝒇 𝑢𝑘 + 𝑩𝒇,𝒅 𝑑𝑘 + 𝑳(𝑦1𝑚𝑒𝑎𝑠,𝑘 − 𝑦̂𝑓,𝑘 ) 𝒙

(32.3)

Jo

The output of the extended model observer 𝑦̂𝑓,𝑘 represents the filtered version of the initial observer output 𝑦̂1,𝑘 . This way the state observer has 𝑛𝑠𝑦𝑠𝑡 + 𝑛𝑓𝑖𝑙𝑡 eigenvalues, where 𝑛𝑠𝑦𝑠𝑡 stands for the system order of initial blender model and 𝑛𝑓𝑖𝑙𝑡 for the order of implemented filter. First 𝑛𝑠𝑦𝑠𝑡 eigenvalues are empirically selected via the observer tests in simulation and at the real plant, and the 𝑛𝑓𝑖𝑙𝑡 ‘filter’ poles are set to zero. Empirically selected poles of the observer are [0.98, 0.98, 0.5]. In order to analyse the observer performance, two observer tests are simulated. In the first observer test, the blender speed remains constant and the inlet mass flow is varied around the nominal value. At t=0s the blender speed and the inlet mass flow are set to their nominal values. At t=80s, the inlet mass flow is increased by 1.5kg/h, and at t=230s it is decreased by 3kg/h. At t=385s it is increased by 4.5kg/h, and at t=555s is set back to the nominal value. In the second observer test, the inlet mass flow is constant and the blender speed varies around the nominal value. At t=0s the blender speed and the inlet mass flow are set to their nominal values. At t=65s, the blender speed is increased by 200rpm, and at t=200s is set to nominal value. Next, at t=325s the blender speed is decreased by 200rpm, and at t=445s it is set back to its nominal value. Figure 16 shows the signals of interest during these observer tests. In both cases, a very good agreement between real and estimated system outputs is achieved.

ro of

Figure 16: Assessment of simulated observer action in observer test 1 (left) and in observer test 2 (right).

Jo

ur

na

lP

re

-p

The MPC performance is examined using an artificially-introduced feeder disturbance. Two experiments are performed with different parameter sets shown in Table 4. Figure 17 provides the simulated and real MPC results in the suggested test, showing the signals computed using the MPC with rate constraint only (left) and the signals computed using the MPC with both rate and magnitude constraints (right). The measurement data for the simulated and measured process are in very good agreement. The first set of parameters allows the controller to keep the blender outlet mass flow very close to the reference value, although a change in the blender speed is not as pronounced as it is using the disturbance feed-forward controller. The results for the MPC with the second parameter set are slightly inferior due to the blender speed constraints introduced into the optimization problem. This test run does not indicate any disadvantages of this controller concept, however, only that the controller ensures that the blender speed remains within the chosen limits.

Figure 17: Assessment of controller action using control strategy 3, parameter set 1 (left) and parameter set 2 (right) in simulation and at the real plant.

5. Conclusions and outlook This paper describes the design and application of several controller strategies for a FBU at a continuous pharmaceutical manufacturing plant. The goal was the damping of brief feeder fluctuations by adjusting the blender speed. Standard controller concepts were compared to a model predictive controller. The first concept suggested (a PI controller with an additional feedback loop) failed to keep the outlet mass flow close to the reference value. The reason was the measurement device setup and the material used. This could be addressed by choosing a material with better properties in terms of mass flow measurement, i.e., a shorter delay due to a lower filter order without corrupting the measurement quality. The results for the second concept (a disturbance feed-forward controller) were satisfactory, i.e., the outlet mass flow was kept close to the reference value. The main advantage of this approach is independence from the outlet mass flow measurements, which would be relatively difficult to accomplish on a continuous production line. In summary, this controller has a very good disturbance rejection and only considers the feeder data. Yet, the actuator constraints cannot be addressed, i.e., rapid changes in the blender speed that could cause a change in the mixing regime cannot be prevented. If these drawbacks are not an issue of the intended application, the second concept is a well suited strategy for outlet mass flow control.

re

-p

ro of

The MPC has many advantages over the other controller concepts. It can consider the input, output and state constraints. By incorporating the blender speed value and rate constraints into the optimization problem, the MPC maintains them within the desired boundaries. By this means an appropriate choice of boundaries prevents any unwanted changes in the mixing regime. Of course, this will also affect the outlet mass flow control performance due to the applied actuator constraints. Additionally, the MPC can handle systems with multiple inputs and outputs. Hold-up is considered as a second system output and its deviation from the reference is controlled via the output error weighting matrix. The MPC tuning is relatively simple, with the controller’s parameters intuitively related to the process performance. However, this controller concept requires a high computational effort, which can be disadvantageous for systems with high sampling frequencies. In addition, the outlet mass flow must be measured in order to estimate the states of the blender model, which may not be feasible on a continuous production line. Consequently, one possible next step would be a test of a trivial observer, which estimates the system’s states solely based on the system’s inputs, i.e., without the outlet mass flow measurement within the MPC framework.

ur

na

lP

On common continuous production lines, a tablet press is placed under the FBU and should receive a constant inlet mass flow with a constant concentration. To date, only the blender’s outlet mass flow and holdup behaviour have been analysed without considering the content uniformity and the concentration of the final product. The next step will be to develop an approach combining the outlet mass flow and concentration control, as suggested in (Singh et al. 2014) and (Singh, Ierapetritou, and Ramachandran 2013). In this constellation, in addition to the blender speed, the feeder’s inlet mass flow set points will be actuating signals. In order to investigate content uniformity at the blender outlet, the existing device setup needs to be expanded by at least one feeder.

Jo

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was partially funded by Austrian Research Promotion Agency FFG, program ”Produktion der Zukunft”, project number 858704.

References: “B&R Industrial Automation.” 2019. B&R, Bus Controller. 2019. https://www.br-automation.com. Bhaskar, Aparajith, Fernando N. Barros, and Ravendra Singh. 2017. “Development and Implementation of an Advanced Model Predictive Control System into Continuous Pharmaceutical Tablet Compaction Process.” International Journal of Pharmaceutics. https://doi.org/10.1016/j.ijpharm.2017.10.003. Chen, C.-T. 1998. Linear System Theory and Design. Oxford University Press. “Coperion GmbH.” 2019. Coperion K-Tron. 2019. https://www.coperion.com/en/products-services/processequipment/feeders/twin-screw-feeders. “Evon GmbH.” 2019. Evon GmbH, XAMControl. 2019. https://evon-automation.com/xamcontrol/. Ferreau, Hans Joachim, Christian Kirches, Andreas Potschka, Hans Georg Bock, and Moritz Diehl. 2014. “QpOASES: A Parametric Active-Set Algorithm for Quadratic Programming.” Mathematical Programming Computation. https://doi.org/10.1007/s12532-014-0071-1. G.F. Franklin, J.D. Powell, M.L. Workman. 2002. Digital Control of Dynamic Systems. Addison-Wesley.

ro of

“HBM GmbH.” 2019. HBM, Catch-Scale. 2019. https://www.hbm.com/en/3020/pw22-high-speed-singlepoint-load-cell-for-dynamic-weighing. “Hosokawa GmbH.” 2019. Hosokawa, Modulo-Mix. 2019. https://www.hosokawa-micronbv.com/technologies/mixing-equipment/continuous-mixing-solutions/modulomix-continuous-modularmixer.html.

-p

Lee, Sau L., Thomas F. O’Connor, Xiaochuan Yang, Celia N. Cruz, Sharmista Chatterjee, Rapti D. Madurawe, Christine M.V. Moore, Lawrence X. Yu, and Janet Woodcock. 2015. “Modernizing Pharmaceutical Manufacturing: From Batch to Continuous Production.” Journal of Pharmaceutical Innovation. https://doi.org/10.1007/s12247-015-9215-8.

re

Maciejowski, J.M. 2002. Predictive Control With Constraints. 1st ed. Pearson Education Limited. “MathWorks GmbH.” 2019. MathWorks, System Identification. 2019. www.mathworks.com/products/sysid.html.

lP

Rehrl, Jakob, Julia Kruisz, Stephan Sacher, Johannes Khinast, and Martin Horn. 2016. “Optimized Continuous Pharmaceutical Manufacturing via Model-Predictive Control.” International Journal of Pharmaceutics. https://doi.org/10.1016/j.ijpharm.2016.06.024.

na

Singh, Ravendra. 2018. “Automation of Continuous Pharmaceutical Manufacturing Process.” In Computer Aided Chemical Engineering. https://doi.org/10.1016/B978-0-444-63963-9.00017-8.

ur

Singh, Ravendra, Dana Barrasso, Anwesha Chaudhury, Maitraye Sen, Marianthi Ierapetritou, and Rohit Ramachandran. 2014. “Closed-Loop Feedback Control of a Continuous Pharmaceutical Tablet Manufacturing Process via Wet Granulation.” Journal of Pharmaceutical Innovation. https://doi.org/10.1007/s12247-014-9170-9.

Jo

Singh, Ravendra, Marianthi Ierapetritou, and Rohit Ramachandran. 2013. “System-Wide Hybrid MPC-PID Control of a Continuous Pharmaceutical Tablet Manufacturing Process via Direct Compaction.” European Journal of Pharmaceutics and Biopharmaceutics. https://doi.org/10.1016/j.ejpb.2013.02.019. “Softing GmbH.” 2019. Softing, Profibus. 2019. https.//industrial.softing.com. Su, Qinglin, Yasasvi Bommireddy, Marcial Gonzalez, Gintaras V. Reklaitis, and Zoltan K. Nagy. 2018. “Variation and Risk Analysis in Tablet Press Control for Continuous Manufacturing of Solid Dosage via Direct Compaction.” In Computer Aided Chemical Engineering. https://doi.org/10.1016/B978-0-44464241-7.50108-7. Su, Qinglin, Mariana Moreno, Sudarshan Ganesh, Gintaras V. Reklaitis, and Zoltan K. Nagy. 2018. “Resilience and Risk Analysis of Fault-Tolerant Process Control Design in Continuous Pharmaceutical Manufacturing.” Journal of Loss Prevention in the Process Industries. https://doi.org/10.1016/j.jlp.2018.07.015.

Table 1: Filter type, cut-off frequency 𝑓𝑐𝑢𝑡−𝑜𝑓𝑓 and order 𝑛𝑓𝑖𝑙𝑡 . 𝑓𝑐𝑢𝑡−𝑜𝑓𝑓 𝐻𝑧 0.0578 0.0578

No. Type

32 32

ur

na

lP

re

-p

ro of

low pass linear phase FIR low pass zero-phase forwardbackward

Jo

1 2

𝑛𝑓𝑖𝑙𝑡

Table 2: Coefficients of identified transfer functions. Transfer function

𝑎0

𝑎1

𝑏0

𝑏1

𝑏2

𝑇𝑑 𝑠

0.03261

0.335

0

0.00457

0.003491

0

𝑃𝑑 (𝑠)

0.2182

0.956

0.2088

0

0

0

𝑃𝑑𝑒𝑙𝑎𝑦 (𝑠)

0.0626

0.5836

0

0.008653

0.007297

8

𝑃𝑑,𝑑𝑒𝑙𝑎𝑦 (𝑠)

0.2324

0.8422

0.2224

0

0

6.65

Jo

ur

na

lP

re

-p

ro of

𝑃(𝑠)

Table 3: State space model parameters. 𝑨 0.0309 0.9239 −0.0322

0.0004 0.0159] 0.9993

𝑩𝒅

𝑪

−4.3042𝑒 − 5 [ 0.0001 ] 5.0884𝑒 − 5

−0.0019 [ 9.8023𝑒 − 5] 0.00023

[

−33.6175 −52.8995

2.0317 31.9043

0.0636 ] −134.042

𝑫

𝑫𝒅

0 [ ] 0

0 [ ] 0

Jo

ur

na

lP

re

-p

ro of

0.9503 [ −0.0356 −0.0098

𝑩

Table 4: MPC parameters, constraints of the blender speed (deviation from the nominal value, Section 2) and the blender speed change.

Run No. 1

𝑛𝑝

40

𝑛𝑐

40

𝑅𝑢,𝑖

𝑸𝒊

𝑅∆,𝒊

size (𝑝 𝑥 𝑝) size(𝑚 𝑥 𝑚) size(𝑚 𝑥 𝑚) 165000 0 [1] [1.5] [ ] 0 200

𝑢𝑚𝑖𝑛

𝑢𝑚𝑎𝑥

∆𝑢𝑚𝑖𝑛

∆𝑢𝑚𝑎𝑥

𝑟𝑝𝑚

𝑟𝑝𝑚

𝑟𝑝𝑚

𝑟𝑝𝑚

-200

200

-10

10

165000 0 [1] [1.5] -70 20 -10 10 ] 0 200 Note: MPC information and detailed derivation of matrices 𝑯, 𝒇𝑻 , 𝑾 and 𝒘 can be found in (Maciejowski 2002). 40

[

ur

na

lP

re

-p

ro of

40

Jo

2