Model predictive control of continuous layering granulation in fluidised beds with internal product classification

Journal of Process Control 45 (2016) 65–75

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Model predictive control of continuous layering granulation in fluidised beds with internal product classification Andreas Bück a,∗ , Robert Dürr b , Martin Schmidt a , Evangelos Tsotsas a a b

Thermal Process Engineering/NaWiTec, Otto von Guericke University Magdeburg, Universitätsplatz 2, D-39106 Magdeburg, Germany Automation/Modeling, Otto von Guericke University Magdeburg, Universitätsplatz 2, D-39106 Magdeburg, Germany

a r t i c l e

i n f o

Article history: Received 7 August 2015 Received in revised form 4 July 2016 Accepted 17 July 2016 Keywords: Layering Particle formation Fluidised bed Population balance modelling Model predictive control

a b s t r a c t Continuous fluidised bed layering with internal product separation is known to possess unstable steady-states, yielding sustained nonlinear oscillations in the particle size distribution, the key product parameter. In this work, the well-known framework of constrained linear model predictive control is applied to stabilise corresponding unstable steady-states. Performance of the controllers is discussed with respect to design parameters and process constraints, as is robustness to overspray generation and size of formed overspray particles. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Fluidised bed technology has found many applications in solids processing, for instance in particle formation and drying [1–3]. In particle formation, fluidised bed technology is mainly used for agglomeration, layering and coating of particles. The main difference between these formation processes is that in agglomeration new particles are formed by coalescence of primary particles. The coalescence can be initiated in different ways, for example the spraying of a binding agent, molecular or electrostatic forces, or thermal effects, e.g., glass transition, melting or sintering [4]. In coating and layering, a suspension or solution is sprayed on the fluidised bed, creating a liquid film on the surface of the particles. By evaporation of the liquid new solid layers are formed, resulting in an onion-like growth of the particles. Whereas coating aims at functionalisation of particles, e.g., taste or odour masking, and therefore only thin layers are applied, layering aims at larger thicknesses, i.e., the formation of large particles of a uniform solid material. The aim in all processes is the production of particles with pre-defined properties, for example the particle size distribution which influences many user-properties, e.g., the re-wettability of agglomerates (and thereby the instant behaviour), or the bulk

∗ Corresponding author. Tel.: +49 391 6718319; fax: +49 391 6718265. E-mail address: [email protected] (A. Bück). http://dx.doi.org/10.1016/j.jprocont.2016.07.003 0959-1524/© 2016 Elsevier Ltd. All rights reserved.

density (via the standard deviation of the size distribution with respect to the mean particle size) which is of huge importance in packaging, transportation and storage of the products [5,6]. Considering layering processes only in the following, from an industrial point of view, continuous operation under steady-state conditions is desired as it allows for a constant product mass flow with constant, user-specified properties. Layering granulation in fluidised beds is usually performed in one of two configurations: (1) internal classification, wherein the product removal is based on a pre-defined mean product particle size by adjusting a gas flow in the outlet tube; or (2) external classification, where particles are removed from the apparatus without classification which is done externally using a screen-and-mill cycle. Both variants can, depending on process parameters which are based on product specifications, show different dynamic behaviour: For some combinations, a stable steady-state is obtained, for others the process is unstable, e.g., yielding sustained oscillations in the particle size distribution in the apparatus. In case of the layering process with external classification, this effect was first described by Heinrich et al. [7], and numerically investigated using bifurcation analysis by Radichkov et al. [8]. Recently, using an advanced process description, the influence of functional zones in the apparatus on process stability was investigated by Dreyschultze et al. [9]. Experimental observation of oscillating behaviour was first reported in the patent of Schütte et al. [10] but without detailed information on the experiments. Schmidt et al. [11] were able to

66

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

confirm these observations repeatedly, for a range of operating conditions. For the configuration with internal classification, Vreman et al. [12], based on a limited number of observations at an industrial plant, came up with a model which was able to show oscillatory behaviour. They explained the different dynamic behaviour with the counterplay of nuclei formation and particle growth in the bed, depending on the distance between the spray nozzle and the particle bed. Recently, Schmidt et al. [13] could show experimentally the different dynamic behaviour for a large set of operating conditions, also investigating the influence of heat and mass transfer on process stability, showing significant influence of the drying conditions on nuclei formation. The oscillatory behaviour leads to an oscillating product mass flow which then requires special attention in all post-processing steps, e.g., drying, cooling, packaging, transportation and storage. Furthermore, there is the risk of overloading the fluidised bed or the danger of emptying it completely, as also shown by Schmidt et al. From the point of process operation there is therefore the need of stabilising these unstable processes in order to achieve the main aim, i.e., operation of the plants at steady-state with a constant product mass flow rate with constant product properties, as well as preventing the drift of the process into dangerous process states as mentioned before. Although there are some works on control of agglomeration (granulation) processes, for instance in drums [14,15], feedback control of layering granulation in fluidised beds has so far only received limited attention in the literature: the works of Cotabarren et al. [16], Palis and Kienle [17–19] and Bück et al. [20] are currently the only focussing on this class of processes. However, most of these works refer to the continuous process with external classification, with the exception of Palis et al. [18] and Bück et al. [20] in which H∞ -loop shaping and robust PI controller theory were employed. In these cases, involved mathematical theory had to be used to derive the feedback law and to provide required robustness properties. Although the resulting controllers can be implemented at a plant, the used theory cannot be considered user-friendly with respect to the plant operator, possibly resulting in permanently switchedoff controllers in case of experience of counter-intuitive controller action. In this contribution, the well-studied concept of constrained linear model predictive control (MPC) is applied for the first time for the design of a suitable feedback controller to stabilise openloop unstable operating points. Furthermore, process constraints are explicitly incorporated in the design, and the performance with respect to practically relevant unmeasurable process uncertainties, overspray generation and size of formed seeds, is studied. The manuscript is structured as follows: first, the process under consideration is described in more detail and a dynamic process model, based on population balances, is derived. Then the model predictive controller is designed, followed by a presentation and discussion of simulation results. Conclusions and outlook then form the final part of this work.

2. Process description and modelling In layering processes, a solid containing liquid, e.g. a solution or suspension, is sprayed onto a bed of fluidised seed particles (Fig. 1). Fluidisation, i.e., a fluid-like behaviour of the particles, is achieved by a heated gas flow which enters the fluidised bed from the bottom through a distributor plate. The heat is used to evaporate the liquid from the particle surface, where the droplets have spread after deposition, so that over time a new solid layer is formed.

Fig. 1. Simplified process scheme of fluidised bed layering granulation with internal classification.

Particles having achieved the desired product size  1 are discharged via a classifying tube which is inserted centrally into the distributor plate. This tube is supplied with a separate gas flow that is adjusted such that the stationary sinking velocity of the particles in the tube corresponds to the product particle size  1 . Ideally, all particles having a size smaller than  1 are blown back into the fluidisation chamber for further growth; in real applications, however, also fractions of the particles with  <  1 are discharged and with  >  1 are blown back, i.e., the classifying tube possesses a separation function. New particles can be supplied in two different ways: externally by adding pre-produced seed particles of a given size, or internally by overspray generation. External production requires a separate process and is in most cases used in coating applications. Overspray is made of seed (nuclei) particles that are produced by pre-drying of sprayed droplets before they come into contact with the particle bed, thus forming individual albeit very small particles. The amount of overspray depends on the drying conditions, e.g., the gas-side conditions, and on the bed characteristics, for example the distance between the bed surface and the spray nozzle and the porosity of the fluidised bed which determine the probability of a droplet-particle collision. Instability of continuous operation in the form of self-sustained oscillations is a counter-play between the discharge of particles and the supply of new particles by overspray generation, as shown by Vreman et al. [12] and Schmidt et al. [13]. Modelling of the layering granulation process follows closely the development given in Vreman et al. [12]—but with several, practice-related comments, modifications and extensions. On a macroscopic scale the particle formation process can be conveniently modelled within the population balance framework [21]. Alternatively, CFD-DEM and Monte-Carlo modelling is applied to particle formation processes [22–32]. They allow investigation of effects on micro- and mesocale, fully resolving single particle interaction. Although powerful, especially to obtain information on process kinetics, they are still too expensive in terms of online computation and thus only population balance modelling is considered here.

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

∂n ∂ = − (Gn) − n˙ discharge (t, ) + n˙ nuc (t, ), ∂t ∂

(1)

1

h < hnoz

∞

 k n(t, ) d,

k (t) =

k = 0, 1, . . .,

(2)

min

which can often be related to physical quantities. In the present case, 0 corresponds to the total number of particles in the bed; 1 to the total length of the particles; 2 is proportional to the total surface area, and 3 is proportional to the total volume of the particles in the bed. The growth rate G is determined by the mass of solid that is sprayed into the apparatus and its distribution on the particles. A growth model, which is commonly used in layering granulation, was given by Mörl et al. [35]: ˙ solid (t) 2 (1 − o(t)) M 2 (1 − o(t)) ˚(t) G= = , part 2 (t) 2 (t)

(3)

which assumes spherical particles of mass density part , a sizeindependent distribution of the sprayed solid on the total particle surface (Abed = 2 ), and the build-up of compact layers. In this equation, the overspray generation is denoted by o. As was explained earlier, a portion of the spray does not contribute to layering growth but to formation of individual nuclei by overspray. Vreman et al. correlated this fraction of the spray to the distance between the bed and spray nozzle, introducing a parameter o* which describes the minimum overspray generation and also incorporates implicitly the influence of the drying conditions. Their correlation is shown in Fig. 2 and—due to the lack of a better description available at the moment—is used in the process model:



o(t) = o∗ + max 0, (1 − o∗ )

(hnoz − hbed (t)) hnoz



,

(4)

where the current bed height is calculated from hbed (t) =

Vbed (t) , (1 − εbed )Aapp

Vbed (t) =

0.6

0.4

0.2

where G is the growth rate of the particles, and n˙ discharge and n˙ nuc represent the particle fluxes due to discharge via the classifying tube and nuclei formation, respectively. Furthermore, the notion of a full-order moment k (t) will be used:



h > hnoz

0.8

overspray fraction o [−]

Introducing the number density distribution of particles, n, with respect to the particle size (diameter), . In the following it is assumed that the fluidised bed is well-mixed, i.e., all spatial dependencies are neglected. Furthermore, only particle growth by layering is considered—agglomeration and breakage effects are neglected. This is achieved by a suitable choice of operating conditions via the Stokes criterion (Ennis et al. [33] and Iveson et al. [34]) as agglomeration and breakage are considered undesired in layering applications. Then, a population balance equation for the temporal evolution of the particle size distribution n(t, ) can be formulated:

67

 3 (t). 6

(5)

The flux of discharged particles, n˙ discharge depends on the characteristics of the outlet tube, as explained earlier. The product diameter  1 corresponds to the stationary sinking velocity which can be manipulated externally by the gas flow in the tube. The discharge function is given by a separation function T(), a strictly increasing function with T( min ) = 0 and max T () = 1, which  ∈ [min ,∞)

takes for ideal separation the shape of a Heaviside step-function, with the step at  =  1 , i.e., T() = H( −  1 ), meaning that all particles with a larger diameter than  1 are discharged. The ideal case is considered in the model of Vreman et al.; in practical applications the separation is usually non-ideal, e.g., due to swarm and cluster formation of particles and the gas flow in the tube, requiring

o = o* 0 0

0.2

0.4

0.6

0.8

1

bed height h [m] Fig. 2. Empirical correlation for overspray generation with respect to bed height, introduced by Vreman et al. [12]. The value hnoz denotes the distance of the nozzle from the distributor plate.

the consideration of usually experimentally determined separation functions T(). Assuming a proportional outlet, the particle flux can be written as: n˙ discharge (t, ) = KT ()n(t, ),

(6)

wherein the drain factor K denotes the probability of a particle entering the outlet tube and is related to the ratio of tube diameter to apparatus diameter. The particle flux due to overspray generation, n˙ nuc can be expressed by the overspray fraction o and the size of an individual, spherical nucleus. Assuming that only monodisperse nuclei of size  nuc are generated (as done by Vreman et al.), the flux can be expressed as n˙ nuc (t, ) =

˙ solid 6oM 3 part nuc

ı( − nuc ) =

6 o ˚(t) 3  nuc

ı( − nuc ),

(7)

where ı denotes the Dirac delta-function. The size of the nuclei,  nuc , is influenced by many different process parameters, for example the initial droplet size, the solid content and also possible formation processes inside the droplets, see, e.g., Handscomb et al. [36], Mezhericher et al. [37], and Bück et al. [38]. In general, the normalised nuclei size distribution can have an arbitrary shape, qnuc (); then the flux can be expressed as: n˙ nuc (t, ) =



∞

6 o ˚(t)

min

 3 qnuc () d

qnuc () .

(8)

Initially, a particle size distribution n(t = 0, ) = f() is specified, as well as a boundary condition, denoting no external seeding: (Gn) (t, min ) = 0 .

(9)

This process model describes the evolution of the particle size distribution in the fluidised bed. As shown qualitatively in Figs. 3 and 4, for different product specifications different process dynamics can be observed in the model (using the simplifying assumptions of Vreman et al.) and also experimentally (data from Schmidt et al. [13]). In order to achieve a stable steady-state operation in all cases, feedback control schemes have to be applied.

68

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

−1

Number density function [m ]

x 10 12 13

x 10 Number density function [m−1]

8 6 4 2 0

10 8

x 10

4

6 4 2 Time [s]

1

0.5

0 0

1.5

2

2.5 x 10

3 −3

2 1.5 1 0.5 0 10

ξ in [m]

3

8

(a) Model simulation.

6

4

x 10

2 −3

4 Time [s]

q [mm −1] 3

x 10

1

2 0

ξ in [m]

0

6 3

(a) Model simulation.

0

16 12

4 0 0

0.2

sep 0.4

0.6

0.8

1

6 3

3

process time t [h]

x

q [mm−1 ]

8

particle size x [mm]

(b) Experimental result.

0

16 12

Fig. 3. Results showing a stable (stationary) normalised volume density distribution in the process chamber (exp. results from Schmidt et al. [13]).

8 3. Controller design

process The population balance model presented in the last section is non-linear and infinite-dimensional from a system-theoretic point of view. For this class of models no standard formalism for feedback controller design is available yet, requiring the use of simplified models derived from the population balance model to calculate appropriate feedback laws. The primary task of the feedback controller is the stabilisation of unstable operating points, especially preventing sustained oscillations, followed by rejection of disturbances. This is to be achieved by fulfilling constraints in the manipulated variables and also showing robustness of the controller with respect to model or parameter uncertainties. Given the product specifications, i.e., the material properties, spraying rate (throughput) and product diameter, the steady-state particle size distribution can be calculated semi-analytically as shown in Appendix A. The steady-state distribution can then be used to linearise the population balance model, resulting in a linear partial differential equation, which in a further step can be discretised with respect to the property coordinate, e.g., by a finite volume method [39], giving a linear approximation of the process dynamics in the vicinity of the steady state: ˜ ˜ ˙ = Ax(t) + Bu(t), x(t)

(10)

where x denotes the state vector consisting of the discretised number density distribution and u is the vector of manipulated variables, expressed as deviations from their corresponding ˜ and B˜ are the usual state space matrices steady-state values. A of appropriate dimensions. Furthermore, some measurements are ˜ available: y = Cx.

x 4 0

time t [h]

0.4

0.2

0

sep 0.6

0.8

1

particle size x [mm]

(b) Experimental result. Fig. 4. Results showing sustained oscillations in the normalised volume density distribution in the process chamber (exp. results from Schmidt et al. [13]).

For practical application of the feedback controllers, a timediscrete formulation is advantageous: x(k + 1) = Ax(k) + Bu(k),

y(k) = Cx(k),

(11)

where k denotes the time instant, as a integer multiple of the sampling time Tsample . It can be written in incremental form, using the increments  · (k) = · (k) − · (k − 1) for u, x and y: x(k + 1) = Ax(k) + BU(k), y(k + 1) = CAx(k) + CBU(k).

(12)

Solving the output equation for y(k + 1) and rearranging the terms an extended model can be derived:



x



y y(k) =



(k + 1) =



0 I

A

0

CA

I

  x y



(k).

x y





(k) +

B CB



U(k)

(13)

(14)

Using the incremental model, it can be seen that integrators for the controlled (and measured) outputs are present in the model

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

formulation. Introducing a new state zT = [(x)T , yT ] the augmented state model can be written in standard form: z(k + 1) = Ae x(k) + Be U(k) ,

y(k) = Ce z(k) ,

(15)

where the matrices are obtained by simple substitution. Choosing the prediction and control horizon as integer multiples of the sampling time Tsample , i.e., Np and Nc , respectively, the output sequence over one prediction horizon can be written more concisely as Y = Fz(k|k) + U ,

(16)

where the vectors Y and U as well as the matrices F and are created by stacking the equations for all sampling times, i.e.,

⎡ ⎢

y(k + 1|k) .. .

Y =⎣

⎤

⎡

⎥ ⎦

U = ⎣

⎢

y(k + Np |k)

⎡

Ce Ae

⎢ ⎢ Ce A2e ⎢ F =⎢ ⎢ .. ⎣ .

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

⎤

U(k)

⎥ ⎦,

.. .

(17)

U(k + Nc − 1)

⎡

Ce Be

⎢ ⎢ Ce Ae Be ⎢ =⎢ .. ⎢ . ⎣

N

Ce Ae p

N −1

Ce Ae p

0

0

Ce Be

0

.. .

Be

N −2

Ce Ae p

.. . Be

⎤ ⎥ ⎥ ⎥ ⎥. ⎥ ⎦

.. . N −Nc

Ce Ae p

Be (18)

The important point is that the sequence can be calculated solely based on the knowledge of the state z(k|k) and the input sequence U(k), . . . U(k + Nc − 1), given U(k) = 0 for k > Nc . This model representation will be used in the following for the design of a model predictive feedback controller to achieve the control aims stated at the beginning of this section. We note that the that the prediction equation can in principle can become numerically ill-conditioned, as for example discussed in Rossiter et al. [40], due to ill-conditioning of the A-matrix and the repeated calculation of matrix powers up to power Np in the calculation of the matrices F and . For the general case, special methods are available to counter this effect, see for example the work of Rossiter et al. [40] and Kouvaritakis and Cannon [41]. After studying the properties of the A-matrix, especially the spread of the absolute values of the largest and smallest eigenvalue, for multiple steady-states of practical interest, revealing only a small spread, an alternative approach is used to reduce numerical illconditioning due to amplification in repeated matrix multiplication which worked sufficiently in all considered cases: the calculation of the powers of A is performed using the ‘binary powering algorithm’ ([42], Ch. 9.2.5) which allows calculation of the kth power of A in at most 2 floor(log 2 (k)) matrix multiplications, for example for Np = 50 at most ten matrix multiplications are required. 3.1. Derivation of model predictive controller The constrained model predictive controller is designed following the standard approaches described in many textbooks, for example [43,41]. As cost functional the quadratic criterion T

J(U) = (R − Y )T (R − Y ) + (U) W (U),

(19)

where W is an input weighting matrix, and R a scaling vector for the reference step trajectory over the prediction horizon, i.e. r = R¯r , where r¯ is a unit-step signal is used. If no further constraints apart from the dynamic state equation are present then an analytic solution for the optimum input sequence (U)opt for each prediction horizon can be obtained:



(U)opt = T + W

−1



TR − T + W

−1

T Fz(k|k).

(20)

69

Here the second part of the right-hand side of the equation can be interpreted as state feedback, the first part is a pre-filter that will guarantee a zero steady-state control error. This becomes more explicit if only (U)opt (k), the actually implemented input, is considered. It can be written as (U)opt (k) = Ky r¯ − KMPC z(k|k),

(21)

and reveals the classical structure of a state feedback controller with a pre-filter for reference tracking. The matrices Ky and KMPC can be obtained from the complete solution over the prediction horizon by taking the first row of the matrices ( T + W ) −1

−1

TR

+ W) respectively. and It can be shown that the optimal solution to the linear unconstrained model predictive control problem is equivalent to the LQ-optimal solution over the same finite time-horizon [44]. However, special care has to be taken to guarantee the stability of the closed-loop system and will be discussed later. Constraints are incorporated into the optimisation problem as equality or inequality constraints; the resulting optimisation problem is a quadratic programme. The quadratic programme cannot be solved analytically and therefore has to be solved on-line by iterative optimisation algorithms, for instance conjugate gradient method, interior point method, or the active-set method are available [45,46]. Input constraints of the form umin ≤ u(k) ≤ umax can be expressed in terms of U by observing that for a single time instant ( T

T F,

u(k) = u(k − 1) + U(k) = Iu(k − 1) + IU(k).

(22)

This yields for the whole control horizon, the following: set of linear inequalities



−(C1 u(k − 1) + C2 U)



C1 u(k − 1) + C2 U



≤

−Umin (k)



.

(23)

Umax (k)

Due to the occurrence of u(k − 1), the constraint matrices have to be updated at every sampling time k and therefore the optimisation problem requires an on-line solution. The complete optimisation problem, a quadratic programme, can thus be stated as T

T

minJ(U) = min (U) ( T +W )(U) − 2(U) T (R − Fz(k|k)), U

U

(24)

subject to

MU ≤ N(k),

(25)

where the constraint matrices have been rearranged into the matrices M and N. The solution to this quadratic programme gives the required control law. Due to the constraints the resulting controller is non-linear, i.e. the closed-loop system is a non-linear dynamic system. 3.2. Closed-loop stability In the unconstrained linear case, closed loop stability of the extended discrete-time state-space model can be investigated by checking the position of the eigenvalues of the closed loop matrix Acl = Ae − Be Kmpc . If all eigenvalue lie in the unit circle in the complex plane, the linear closed-loop system is (asymptotically) stable. This also means that the nonlinear system is (asymptotically) stable as long as the linear approximation represents the dynamics of the nonlinear plant [47,44]. In the constrained case, the situation is more complex as the controller is now, in general, non-linear, yielding a non-linear closed-loop system. Consequently, to investigate closed-loop stability nonlinear methods have to be applied. Two main approaches

70

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

can be found in the literature, based on the optimisation programme for the cost functional J(U) =

Np 

T

T

[z(k) Qz(k) + (U(k)) W (U(k))],

0.1 0.08 0.06

(26)

0.04

k=1

(27)

y(k) = Ce z(k),

(28)

with additional appropriately formulated constraints, for instance on the control increments U. One approach draws from the solution of the unconstrained case which is equivalent to linear-quadratic (LQ) optimal control: stability of the closed-loop system in the presence of constraints can be obtained for sufficiently large prediction horizon (Np → ∞) provided that the pair of matrices (Ae , Q1/2 ) is observable for positive definite matrices Q. In this case, Np can be used as a tuning parameter to achieve closed-loop stability. The other approach uses a terminal weight in the finite-time cost functional to account for the cost on the infinite-time horizon. This idea also draws from the fact that unconstrained linear MPC and LQ-optimal control are equivalent. The terminal weight Q¯ is introduced as



Np −1

J(U) =

T

T



z(k) Qz(k) + (U(k)) W (U(k)) + z(Np )T Q¯ z(Np ).

k=1

(29) For the actual computation of Q¯ a Lyapunov equation can be derived, details are given for instance in Mayne et al. [44]. Closedloop stability can then be shown using the optimal cost functional as a Lyapunov function candidate. Tuning the prediction horizon (Np ), the constrained system is steered into an invariant region of the state space where all constraints are and remain satisfied. From then on the system is unconstrained and stability is guaranteed by the terminal weight that corresponds to an LQ-optimal controller [47,44]. 4. Results and discussion The process model presented was implemented in Matlab on a standard desktop PC, using a discretisation of the partial differential equation by a first-order (upwind) finite volume scheme. The quadratic optimisation programme was solved in each time interval by Matlab’s quadprog routine using an active set strategy. The nominal process conditions for which the controllers were tested, are given in Table A.1 and are the values given in the work of Vreman et al. [12], corresponding to unstable steady-point exhibiting sustained oscillations in the number density distribution (stable limit cycle). Although results are only shown for the ‘Vreman scenario’, similar results can be obtained also for arbitrary separation functions and nuclei size distributions, as the required steady-state distributions can always be calculated as shown in Appendix A. As manipulated variable, the solid volume spray rate ˚ is used as it has direct influence on the stability of the process. Furthermore, it can be easily adjusted (within bounds) at operating plants. Input constraints have to be taken into account, to stay in the coating regime and not to change other particle characteristics, e.g., layer porosity [48]. As measured and controlled quantity the total surface area of particles in the bed Abed =  2 is used. The value can be obtained by simultaneously measuring the normalised particle size distribution and the total bed mass mbed in the apparatus. The bed mass can be obtained easily from pressure drop measurements; for the measurement of normalised particle size distributions several online devices are available, e.g., the IPP-70S (Parsum GmbH, Germany).

0.02 Im λ

z(k + 1) = Ae z(k) + Be U,

0 −0.02 −0.04 −0.06 −0.08 −0.1 0.8

0.85

0.9

Re λ

0.95

1

1.05

Fig. 5. Distribution of eigenvalues of the extended process model at the steadystate corresponding to process parameters in Table A.1. Crosses denote eigenvalues, the dashed line is part of the unit circle, marking the stability boundary (stable: left of the boundary).

In all cases, the cost functional is given as T

J(U) = (R − Y )T (R − Y ) + (U) W (U),

(30)

where the weighting matrix W is diagonal. A first parameterisation is done empirically (‘Bryson’s rule’), however, the influence of this weight on the performance is also investigated in the following. 4.1. Stabilisation of unstable steady-states First, the unstable operating point resulting from the process conditions given in Table A.1 is to be stabilised. The eigenvalue distribution of the extended linearised open-loop process model is shown in Fig. 5. A pair of complex eigenvalues outside the unit circle is observed, resulting in an unstable steady-state. The manipulated variable is constrained to ±20% of the steadystate value ˚. The diagonal elements of the matrix W are set initially to 108 and are later adjusted to 1012 . The initial displacement of the process is set to 2% of the steady-state distribution, making sure that the linear model still approximates the nonlinear process. Closed-loop stability is realised by the terminal constraint in the problem formulation as described earlier. The length of the prediction horizon Np is based on the eigenvalue closest to the stability boundary, but is one of the tuning variables in the MPC formulation. The control horizon Nc is set to Nc = Np /3, rounded to the nearest integer. Sampling times investigated are Tsample = {60, 300} s, the corresponding lengths of the prediction horizon are Np = 121 and Np = 25, respectively. Fig. 6 shows the evaluation of the measured and controlled variable y(kTsample ) = 2 (kTsample ) over process time under application of the model predictive controller for the two sampling times and the two control weights. The labelling of the curves is ‘sampling time/log 10 weight’, e.g., ‘60/8 for a sampling time of 60s and diagonal weighting matrix W with all diagonal elements set to 108 . It can be seen that for all tested controller configurations, the steady-state value is achieved. The performance (speed) depends on the sampling time and the weighting: while the combination ‘60/8’ is faster than ‘60/12’—due to the fast interaction with the process and the larger region of allowed input values—the combination ‘300/8’ achieves the result later than the configuration ‘300/12’. The reason for this is that the lower weighting allows for larger input signals which, when applied for too long without

A. Bück et al. / Journal of Process Control 45 (2016) 65–75 12

2780

2720 2700 2680 2660

Reference uncontrolled 60/8 60/12 300/8 300/12

10 −1

2740

x 10

Number density function [m ]

steady−state 60/8 60/12 300/8 300/12

2760

Controlled variable [m2]

71

8

6

4

2 2640 0

2

4

6 Time [h]

8

10

12

Fig. 6. Measured and controlled variable y(kTsample ) = 2 (kTsample ) for different sampling times and control weights. The labelling of the curves is ‘sampling time/log 10 weight’.

0

0.5

1

1.5 ξ in [m]

2

2.5 −3

x 10

Fig. 8. Achieved number density distributions close to the end of the simulation interval (approx. 11.5 h) for the different combinations in sampling time and control weights.

−4

x 10

18 steady−state 60/8 60/12 300/8 300/12

3

−1

Manipulated variable [mm s ]

1.6 1.5 1.4 1.3

16 Normalised error E2

1.7

14 12

uncontrolled 60/8 60/12 300/8 300/12

10 8 6 4 2

1.2

0 0

1.1

0

2

4

6 Time [h]

8

10

12

Fig. 7. Input profiles generated by the model predictive controller for different sampling times and control weights.

update (e.g., 300 s), drive the process temporarily farther away from the set-point, as can be observed in the larger maximum amplitude in Fig. 6. Fig. 7 shows the corresponding input profiles for the different combinations of sampling time and control weight. Here, for the combination ‘300/8’ a high variability in the input signal can be observed, initially oscillating strongly between the input constraints and resulting in the slower performance of this controller parameterisation. The results for the number density distribution for the different combinations are shown in Fig. 8 close to the end of the simulated time interval. The steady-state profile is attained in all scenarios. In order to evaluate the performance of the controller, Fig. 9 shows the (discretised) L2 -error norm of the number density distribution with respect to the steady-state distribution normalised to the initial value. A similar performance in the state error can be observed when compared to the performance in controlled variable. In all cases, the error goes to zero, i.e., the steady-state distribution is reached.

2

4

6 Time [h]

8

10

12

Fig. 9. Plot of the error in number density distribution with respect to the steadystate distribution measured in the L2 -norm, normalised to the initial value of the error.

In all cases, a stabilisation of the unstable steady-state for initial state deviations can be achieved. For increased sampling time, the input weights have to be increased to keep the process in the vicinity of the steady-state without decrease of controller performance. The computational time required for the solution of the quadratic optimisation programme was small compared to the sampling time in every case (maximum: <1%), i.e., a faster-than-realtime implementation at the process plant is possible. A comparison of the designed model predictive controller with a standard LQ-optimal controller using the same state-space model and weighting was performed. In the unconstrained case, i.e. no bounds on the manipulated variable, the controller generated physically irrelevant values, for example non-positive spray rates. In order to prevent this, the bounds on the manipulated variable were implemented by clipping the values at the minimum and maximum, respectively. In order to achieve closed-loop stability, the input weighting had to be increased significantly, still resulting in stronger oscillatory behaviour in manipulated variable when compared with the profiles generated by the model predictive controller. This is shown in Figs. 10 and 11 where exemplary results are presented for a sampling of 60 s and an input weighting of 1015 . Although a stabilisation is achieved, i.e. the controlled variable achieves the unstable steady-state value, significant controller

72

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

35

2750 Reference controlled

2730 Controlled variable [m2]

uncontrolled 60/12 300/12

30 Normalised error E2

2740

2720 2710 2700 2690

25 20 15 10 5

2680 0 0

2670 2660 2650 0

2

4

6 t in [h]

8

10

12

2

4

6 Time [h]

8

10

12

Fig. 12. Plot of the error in number density distribution with respect to the steadystate distribution measured in the L2 -norm for the case of a plant mismatch of +10% in the nuclei size  nuc formed by overspray.

14

Fig. 10. Plot the controlled variable under LQ-optimal control using a sampling time of 60 s and an input weighting of 1015 .

12 −4

Reference Controller

3

Manipulated variable [mm s−1]

1.4

1.35

Normalised error E2

x 10

10 8 6 4

1.3 uncontrolled 60/12 300/12

2 1.25 0 0

5

10

15

20 Time [h]

25

30

35

1.2

0

2

4

6 t in [h]

8

10

12

Fig. 11. Plot the manipulated variable under LQ-optimal control using a sampling time of 60 s and an input weighting of 1015 . Bounds on the input are realised by clipping; in this example bounds are not active due to the high input weighting.

action is required. The resulting dynamic closed-loop behaviour is comparable to the model predictive case ‘300/12’, however, there less control action is required. For that reason, in the following only the constrained model predictive controller is considered. 4.2. Model mismatch The main assumption in the previous investigations is that the process model is perfect, however, this cannot be guaranteed in the actual application. Therefore, the performance of the model predictive controller is evaluated for model mismatch in the vicinity of the unstable operating point. For this study, in addition to the initial displacement from the steady-state by 2%, mismatch in the size of seed particles generated by overspray,  nuc , and in the minimum value of overspray generation o* are considered. Both are closely related to external conditions, e.g., heat and mass transfer between the particulate and the fluid phase. In the study, for both variables deviations of +10% are considered in the nonlinear model, i.e., the linearised model

Fig. 13. Plot of the error in number density distribution with respect to the steadystate distribution measured in the L2 -norm for the case of a plant mismatch of +10% in the steady-state overspray formation rate o*.

used for controller design is the same as in the pure stabilisation study. As can be seen, in Fig. 12 for mismatch in  nuc , Fig. 13 for mismatch in o*, and Fig. 14 for both simultaneously, even in the presence of model mismatch, the unstable steady-state can be stabilised by the model-predictive controller. This is again shown by the error in number density distribution with respect to the (mismatch-free) steady-state distribution which tends in all cases to a non-zero but steady finite value, whereas in the uncontrolled case the error shows oscillatory behaviour, i.e., oscillations in the number density distribution occur. The occurrence of a non-zero steady-state error is not due to a shortcoming of the MPC algorithm as such, but due to the available manipulated variable u = ˚ and the physics of the process. Consider first the uncertainty in the size of the formed seed (nuclei): from the analytic solution given in the Appendix, it can be seen that the steady-state distribution is parameterised by the sizes  nuc ,  1 and the separation efficiency. The latter one is the same in all cases, i.e. the steady-state distribution is pinned between  nuc and  1 , as the formed seeds are the particles having smallest size in the cases under consideration. The manipulated variable, the spray rate, can only influence the magnitude of the distribution between these two sizes, i.e. it cannot compensate for the change in seed size. Therefore, the process tends to a new stabilised steady-state (with

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

73

−4

30 1.65 25

x 10

20

3

Normalised error E2

−1

Manipulated variable [mm s ]

1.6

15

10

5 uncontrolled 300/12 0 0

5

10

15

20 Time [h]

25

30

35

Fig. 14. Plot of the error in number density distribution with respect to the steadystate distribution measured in the L2 -norm, for plant mismatch of +10% in both,  nuc and o*.

a new steady-state distribution). The difference between the two steady-state distributions (measured as an integral value) attains a constant non-zero value. More explicitly: the original steady-state is not reachable under this combination of manipulated variable and uncertainty. Although the error does not vanish, it is bounded, achieving a stabilisation of the nonlinear process, in contrast to open-loop operation which shows sustained nonlinear oscillations. The actual difference is also influenced by the weighting in the cost functional, i.e. with very low weighting on the manipulated variable and very strict weighting on the state, the L2 -error measure could be reduced (by increasing u) until constraints on the manipulated variable become active. However, even in that case the achieved steady-state distribution would differ as a plot of the L∞ -norm of the difference would reveal. In terms of the minimum overspray generation o*, the ability of the controller to cope with the uncertainty is limited by the diametral interaction of the manipulated variable u = ˚ with a) seed formation (interaction o˚) and b) the overall particle growth rate G (interaction (1 − o)˚). In order to compensate an increase in the minimum overspray generation, the spray rate has to be decreased. This, however, yields a decrease in overall growth rate G which in turn increases the time required for particles to grow from  nuc to  1 , thereby increasing the hold-up in the apparatus. This increased hold-up then results in an increased magnitude of the number density distribution in the size range from  nuc to  1 with respect to the nominal steady-state number density distribution, leading to the observed offset in the L2 measure of the error in the controlled distribution. Vice versa, in order to minimise this error, the spray rate is increase until the upper constraint becomes active, as actually happens in the present case shown for the manipulated variable in Fig. 15. Whereas model mismatch in  nuc is adjusted for quite fast independent of the tested sampling times, it takes much longer for a model mismatch in o* to be adjusted, again highlighting the sensitivity of the process with respect to overspray generation. However, in case of simultaneous mismatch (Fig. 14), the time behaviour is again dominated by the model mismatch in  nuc . In all scenarios, the process is steered to a new steady-state, which is nevertheless stabilised if open-loop unstable. This requires of course that the new steady-state is still in the region of validity of the linear approximation of the nonlinear process dynamics and input constraints remain inactive.

1.55 1.5

steady−state 60/12 300/12

1.45 1.4 1.35 1.3 0

2

4

6 Time [h]

8

10

12

Fig. 15. Controller output in case of unmeasured model mismatch in minimum overspray generation o*. (Explanation of this behaviour is given in the text.)

The computational times required for the solution of the optimisation programme are again less than 1% of the sampling time, allowing real-time implementation of the calculated control inputs at process plants. 5. Conclusions and outlook In this work, it was shown that constrained linear model predictive control is a viable option for stabilisation of open-loop unstable operating points in continuous fluidised bed layering with internal product separation. It was further shown that the controller can be designed directly from product specifications using a semianalytical expression for the required steady-state particle size distribution. The performance of the controller depends on both, input weighting and sampling time, using a terminal constraint, however, closed-loop stability can be guaranteed. Finally, it was shown that the model predictive controller possesses a certain robustness with respect to model mismatch, steering the process to new stabilised steady-states. Further extensions of this concept can be: design of a start-up strategy based on a set of operating points and controllers by gain scheduling, or extension to full nonlinear model predictive control covering all process conditions utilising the result of the linear model predictive controller as initial solution. From the process point of view, the influence of agglomeration as seed generating step should be investigated as it was pointed out to be important for process performance [13], also the formation of functional zones can be considered, as well as the inclusion of other processes influencing the layering process, for example heat and mass transfer. Acknowledgements The authors gratefully acknowledge the funding of this work by the German Federal Ministry of Science and Education (BMBF) as part of the InnoProfile-Transfer project NaWiTec (03IPT701X). Appendix A. Derivation of steady-state particle size distribution In the following the steady-state particle size distribution required for controller design is derived for the general process model given by Eqs. (1–9). For the case of ideal separation in the outlet tube and formation of monodisperse nuclei, this was first

74

A. Bück et al. / Journal of Process Control 45 (2016) 65–75

Table A.1 Process parameters used in simulation studies as reported by Vreman et al. [12]. Aapp K ˚0 ˚ hnoz o* εbed  nuc 1

5.0 1.92 × 10−4 1.67 × 10−4 0.8 ˚0 0.44 0.028 0.55 0.3 × 10−3 0.7 × 10−3

m2 s−1 mm3 s−1 mm3 s−1 m – – m m

dns () − 0 = −Gs () d





dGs () + KT () d

ns () + n˙ nuc,s ()

(A.6)

the steady-state number density distribution can be obtained by first by calculating a homogeneous solution, followed by calculation of a particular solution, e.g. by variation of constants:

Table A.2 Main notation. hbed hnoz n o q3 u Aapp B0 A, B, C, F G J K Nc Np Q R Tsample T W Y εbed

j  1  nuc ˚

where B0,s accounts for possible external steady-state supply of monodisperse nuclei of size  min . Note, that the Vreman model is recovered by setting 3 )) ı( −  and

2 = 1, B0,s = 0, n˙ nuc,s () = ((6 os ˚s )/( nuc nuc ), T() = H( −  1 ). From the steady-state balance

Bed height (m) Distance nozzle-distributor plate (m) Number density function (m−1 ) Mass fraction of solid spray that creates nuclei Normalised number density function (m−1 ) Manipulated variable Apparatus cross-sectional area (m2 ) Number flow rate of external seeds (s−1 ) State space matrices Particle growth rate (m s−1 ) Cost functional Drain factor (s−1 ) Number of control steps Number of prediction steps Weighting matrix Reference over prediction horizon (m2 ) MPC sampling time (s) Separation function Controller design matrices Measurements over prediction horizon (m2 ) Bed porosity Weight jth total moment of number density function (mj ) Particle diameter (m) Product particle diameter (m) Nuclei particle diameter (m) State space matrix Volume spray rate (mm3 s−1 )

ns () =

B0,s exp(−Q ()) + exp(−Q ()) Gs (min )









Bs (z) exp(Q (z)) dz, min

(A.7)

where As () = dGs ()/d + KT () /Gs (), Bs () = n˙ nuc,s ()/Gs () and Q () =



min

As () d.

Note, that due to the occurrence of ns () in the definition of Gs (), the equation above determines the shape of ns only; quantitative information can, however, be obtained by solving simultaneously the equation 0 = Gs () −

K 

k Gk,s (ns ())

(A.8)

k=0

by iteration, for example using a Newton-type method. The shape of the steady-state number density distribution has to be evaluated numerically in most cases due to the possible complexity of Gs () and T(), however, in case of the a propertyindependent growth rate and ideal separation T() = H( −  1 ), the shape of the steady-state distribution can also be evaluated explicitly, as shown by Vreman et al. [12]. References

done by Vreman et al. [12]. We extend the calculation to a larger set of process models considering arbitrary separation functions and nuclei size distributions as well as more general growth laws of the form [48] (Tables A.1 and A.2): G(t, ) =

K 

k Gk (t, ).

(A.1)

k=0

For mass conservation, the restriction K 

k = 1

(A.2)

k=0

is imposed; the growth rates Gk are given by Gk (t, ) =

2 (1 − o)˚(t) k−2 .   k (t)

(A.3)

This extended growth model can be utilised to describe dynamic effects, i.e., dispersion of the particle size distribution (widening, narrowing), which cannot be described by the growth rate introduced in the main text, which corresponds to the case 2 = 1. The steady-state population balance equation then reads: 0=−

d (Gs ()ns ()) − KT ()ns () + n˙ nuc,s (), d

(A.4)

with the boundary condition (Gs ns ) (min ) = B0,s ,

(A.5)

[1] L. Reh, Trends in research and industrial application of fluidization, Verfahrenstechnik 11 (1977) 381–384. [2] D. Kunii, O. Levenspiel, Fluidization Engineering, Butterworth-Heinemann, Boston, 1991. [3] J. Burgschweiger, E. Tsotsas, Experimental investigation and modelling of continuous fluidized bed drying under steady-state and dynamic conditions, Chem. Eng. Sci. 57 (2002) 5021–5038. [4] A. Bück, E. Tsotsas, K. Sommer, Size enlargement, in: B. Elvers (Ed.), Ullmann’s Encyclopedia of Industrial Chemistry, Wiley-VCH, Weinheim, 2014, pp. 1–77. [5] H. Rumpf, A. Gupte, Einflüsse der Porosität und Korngrößenverteilung im Widerstandsgesetz der Porenströmung, Chem. Ing. Tech. 43 (6) (1971) 367–375. [6] E. Tsotsas, Heat and mass transfer in packed beds with fluid flow (M7), in: VDI Heat Atlas, Springer-Verlag, Berlin, 2010, pp. 1327–1341. [7] S. Heinrich, M. Peglow, M. Ihlow, L. Mörl, Particle population modeling in fluidized bed spray granulation – analysis of the steady state and unsteady behavior, Powder Technol. 130 (2003) 154–161. [8] R. Radichkov, T. Müller, A. Kienle, S. Heinrich, M. Peglow, L. Mörl, A numerical bifurcation analysis of continuous fluidized bed spray granulation with external product classification, Chem. Eng. Process. 45 (2006) 826–837. [9] C. Dreyschultze, C. Neugebauer, S. Palis, A. Bück, E. Tsotsas, S. Heinrich, A. Kienle, Influence of zone formation on stability of continuous fluidized bed layering granulation with external product classification, Particuology 23 (2015) 1–7. [10] R. Schütte, A. Ruhs, I. Pelgrims, C.-J. Klasen, L. Kaiser, Verfahren zur Herstellung von Granulat mit periodisch oszillierender Korngrößenverteilung und Vorrichtung zu seiner Durchführung, German Patent DE 19639579 C1, 1963. [11] M. Schmidt, C. Rieck, A. Bück, E. Tsotsas, Experimental investigation of process stability of continuous spray fluidized bed layering with external product separation, Chem. Eng. Sci. 137 (2015) 466–475. [12] A. Vreman, C. van Lare, M. Hounslow, A basic population balance model for fluid bed spray granulation, Chem. Eng. Sci. 64 (2009) 4389–4398. [13] M. Schmidt, A. Bück, E. Tsotsas, Experimental investigation of process stability of continuous spray fluidized bed layering with internal separation, Chem. Eng. Sci. 126 (2015) 55–66.

A. Bück et al. / Journal of Process Control 45 (2016) 65–75 [14] M. Pottmann, B. Ogunnaike, A. Adetayo, B. Ennis, Model-based control of a granulation system, Powder Technol. 108 (2000) 192–201. [15] T. Glaser, C. Sanders, F. Wang, I. Cameron, J. Litster, J.-H. Poon, R. Ramachandran, C. Immanuel, F. Doyle III, Model predictive control of continuous drum granulation, J. Process Control 19 (2009) 615–622. [16] I. Cotabarren, D. Bertin, V. Bucala, J. Pina, Feedback control strategies for continuous industrial fluidized-bed granulation process, Powder Technol. 283 (2015) 415–432. [17] S. Palis, A. Kienle, Stabilization of continuous fluidized bed spray granulation with external product classification, Chem. Eng. Sci. 70 (5) (2012) 200–209. [18] S. Palis, A. Kienle, H∞ loop shaping control for continuous fluidized bed spray granulation with internal product classification, Ind. Eng. Chem. Res. 52 (2012) 408–420. [19] S. Palis, A. Kienle, Discrepancy based control of particulate processes, J. Process Control 24 (2014) 33–46. [20] A. Bück, S. Palis, E. Tsotsas, Model-based control of particle properties in fluidised bed spray granulation, Powder Technol. 270 (2015) 575–583. [21] D. Ramkrishna, Population Balances: Theory and Application to Particulate Systems in Engineering, Academic Press, New York, 2000. [22] A. Di Renzo, F. Di Maio, R. Girimonte, B. Formisani, DEM simulation of the mixing equilibrium in fluidized beds of two solids differing in density, Powder Technol. 184 (2008) 214–223. [23] C. Radeke, B. Glassner, J. Khinast, Large-scale powder mixer simulations using massively parallel GPU architectures, Chem. Eng. Sci. 65 (2010) 6435–6442. [24] L. Fries, S. Antonyuk, S. Heinrich, S. Palzer, DEM-CFD modelling of a fluidized bed spray granulator, Chem. Eng. Sci. 66 (2011) 2340–2355. [25] M. Börner, M. Peglow, E. Tsotsas, Derivation of parameters for a two compartments population balance model of Wurster fluidised bed granulation, Powder Technol. 238 (2013) 122–131. [26] Z. Li, J. Kessel, G. Grünewald, M. Kind, Coupled CFD-PBE simulation of nucleation in fluidized bed spray granulation, Drying Technol. 31 (2013) 1888–1896. [27] K. Lee, T. Matsoukas, Simultaneous coagulation and break-up using constant-N Monte Carlo, Powder Technol. 110 (2000) 82–89. [28] A. Maisels, F. Einar Kruis, H. Fissan, Direct simulation Monte Carlo for simultaneous nucleation, coagulation, and surface growth in dispersed systems, Chem. Eng. Sci. 59 (2004) 2231–2239. [29] H. Zhao, A. Maisels, T. Matsoukas, C. Zheng, Analysis of four Monte Carlo methods for the solution of population balances in dispersed systems, Powder Technol. 173 (2007) 38–50. [30] D. Meimaroglou, C. Kiparissides, Monte Carlo simulation for the solution of the bi-variate dynamic population balance equation in batch particulate systems, Chem. Eng. Sci. 62 (2007) 5295–5299.

75

[31] M. Dernedde, M. Peglow, E. Tsotsas, A novel, structure-tracking Monte Carlo algorithm for spray fluidized bed agglomeration, AIChE J. (2011). [32] M. Dernedde, M. Peglow, E. Tsotsas, Stochastic modelling of fluidized bed granulation: influence of droplet pre-drying, Chem. Eng. Technol. 34 (2011) 1172–1176. [33] B. Ennis, G. Tardos, R. Pfeffer, A microlevel-based characterization of granulation phenomena, Powder Technol. 65 (1991) 257–272. [34] S. Iveson, J. Litster, K. Hapgood, B. Ennis, Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review, Powder Technol. 117 (2001) 3–39. [35] L. Mörl, S. Heinrich, M. Peglow, Fluidized bed spray granulation, in: Granulation, Elsevier B. V, Amsterdam, 2007, pp. 21–188. [36] C. Handscomb, M. Kraft, A. Bayly, A new model for the drying of droplets containing suspended solids, Chem. Eng. Sci. 64 (2009) 628–637. [37] M. Mezhericher, A. Levy, I. Borde, Heat and mass transfer of single droplet/wet particle drying, Chem. Eng. Sci. 63 (2008) 12–23. [38] A. Bück, M. Peglow, M. Naumann, E. Tsotsas, Population balance model for drying of droplets containing aggregating nanoparticles, AIChE J. 58 (11) (2012) 3318–3328. [39] R. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, Basel, 1992. [40] J. Rossiter, B. Kouvaritakis, M. Rice, A numerically robust state-space approach to stable-predictive control strategies, Automatica 34 (1) (1998) 65–73. [41] B. Kouvaritakis, M. Cannon, Model Predictive Control: Classical, Robust, Stochastic, Springer International Publishing, 2016. [42] G. Golub, C. Van Loan, Matrix Computations, fourth ed., The Johns Hopkins University Press, 2013. [43] J. Maciejowski, Predictive Control with Constraints, Prentice Hall, 2001. [44] D. Mayne, J. Rawlings, C. Rao, P. Scokaert, Constrained model predictive control: stability and optimality, Automatica 36 (2000) 789–814. [45] P. Gill, L. Jay, M. Leonard, L. Petzold, V. Sharma, An SQP method for the optimal control of large-scale dynamical systems, J. Comput. Appl. Math. 120 (2000) 197–213. [46] J. Nocedal, S. Wright, Numerical Optimization, Springer Science + Business Media, LLC, London, 2006. [47] H. Chen, F. Allgöwer, Quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability, Automatica 34 (10) (1998) 1205–1217. [48] C. Rieck, T. Hoffmann, A. Bück, M. Peglow, E. Tsotsas, Influence of drying conditions on layer porosity in fluidized bed spray granulation, Powder Technol. 272 (2015) 120–131.