A Mathematical Framework for Spatio-Temporal Control in Industrial Drying
Lennon Ó Náraigh






Corresponding author.  Email: onaraigh@maths.ucd.ie
School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland






We introduce two models of industrial drying - a simplified one-equation model, and a detailed three-equation model.  The purpose of the simplified model is rigorous validation of numerical methods for PDE-constrained optimal control.  The purpose of the detailed model is to be able to predict and control the behaviour of an industrial disk drier.  For both models,  we introduce a fully validated numerical method to compute the optimal source term to maintain the outlet temperature as close as possible to the set-point temperature.  By performing simulations using realistic parameters for industrial driers, we illustrate potential applications of the method.







A Mathematical Framework for Spatio-Temporal Control in Industrial Drying
    Lennon Ó Náraigh
    today
=========================================================================







§ INTRODUCTION



In this article, we consider the mathematical modelling of  industrial driers.  Such machines  are used in the food-processing industry to remove moisture from raw materials; they require large quantities of energy to raise the temperature of the moisture and to evaporate the moisture in the raw materials.  The  machine we have in mind for modelling is the disk dryer, however, the mathematical framework we develop can be applied to a wide variety of driers and heaters.  Such machines do not always operate in steady state, as the stream of input raw materials may vary in temperature and composition.  Therefore, in order for the final moisture level and temperature at the outlet of the machine to have a fixed value, a heat source must be present to deliver a controlled amount of thermal energy to the system, to correct for the inlet fluctuations and hence, to `steer' the system to a target state at the outlet.  Mathematical models for industrial drying can give important insights into the correct way to control such systems.  Hence, the aim of this work is to develop such models, with a view to formulating a control theory to maintain the outlet temperature at a specific set-point.  



 §.§ Literature Review


The simplest possible models in the drying literature are lumped-parameter models.  Here, either the spatial variations or the temporal variations are averaged out, leading to a set of simple ordinary differential equations based on mass- and energy-balances<cit.>.    Examples include Reference <cit.> for a disk dryer, and Reference <cit.> for a rotary dryer.  In these cases, the spatial degrees of freedom are averaged out, leading to a system of ordinary differential equations in the time variable.  The authors then use model predictive control to maintain the outlet variables at a set point.


And yet, often, the temperature and moisture profiles along the length of the machine are important – for instance, it is desirable for the outlet temperature to be at a set point but also, for the temperature inside the drier to stay within certain bounds, for energy efficiency but also, for food hygiene.  For this reason, Reference <cit.> looks at a set of coupled ordinary differential equations for a rotary drier, with each set of equations representing a separate segment of the rotary drier.  The authors look at how the temperature profile along the drier depends on inlet conditions and flow rate – effectively the reverse of the previous examples where now the spatial variable (i.e. the distance x from the inlet, as opposed to time) is the key variable of interest.  The authors also look at how the system responds to a step change in the dryer speed.

In the present work, both the temporal and spatial variables are important.  Hence, we look at a partial-differential-equation (PDE) model for industrial driers.  In this way, the hyperbolic nature of such systems becomes apparent – information from the inlet is convected to the outlet.  Hence, controlling the moisture level and the temperature at the outlet requires the propagation of information from the inlet, as well as manipulation of the heat source inside the drier, to `steer' the system towards a desired outlet state.

This introduces the idea of optimal control theory.  Optimal control theory is well understood in the context of ordinary differential equations, where it is required to minimize a cost functional based on state variables, subject to the constraint that the state variables satisfy a given mathematical model.  
A time-dependent forcing term (the control variable) provides the necessary additional degrees of freedom so that the minimum can be realised – the control variable which realises the minimum of the cost functional is called the optimal control.  

For ordinary differential equations, the optimal control can be computed using Pontryagin's Maximum Principle<cit.> – the idea here is to introduce a constrained penalty function comprising the basic objective function to be minimized, as well as a set of constraints (parametrized by Lagrange multipliers) enforcing the condition that the state variables satisfy the mathematical model.  By taking variations of the constrained functional, one obtains a set of forward equations (the  mathematical model itself), as well as a set of adjoint equations (involving the Lagrange multipliers), the solution of which gives the optimal control. 


In the present case, we apply optimal control theory to partial differential equations, where a similar technique applies.
The reader is referred to Reference <cit.> for details.  The same technique has been used in a PDE model for the optimal control of chemical reactors in microchannels<cit.>, a PDE model for the optimal draw speed in the production of glass (`tube drawing')<cit.>, and PDE model for the optimal air temperature / air speed for melt spinning, a key manufacturing technique in the production of synthetic textile fibers<cit.>.  The technique has also been applied to Burgers equation<cit.>, a canonical non-steady problem and similar in its mathematical form to the models considered in the present work.




 §.§ Plan of the paper


References <cit.> and <cit.> use generalized Newton methods for the minimization of the constrained cost functional and the computation of optimal controls, an approach we also use in this paper.  However, our own first steps are more cautious, and we start with a simpified model of industrial heating (Section <ref>).  Crucially, the simplified model admits exact analytical solutions – including for the  optimal control.  
In Section <ref>, we introduce a numerical optimal-control theory for the simplified model, and we use the available  analytical solutions to validate the numerical method rigorously.




Next, we introduce a more detailed model of an industrial drier, which convects the product from the inlet to the outlet (Section <ref>).  The product is made up of two streams: dry solid and liquid.  The model allows for the evaporation of the liquid and also involves the temperature of the streams – hence, a three-equation model.  

We  look at the equilibrium solutions of the model, corresponding to steady inlet conditions, as well as the linear stability of such equilibrium solutions (Section <ref>).  We investigate small-amplitude perturbations away from equilibrium, and demonstrate analytically that these perturbations can be suppressed so that the outlet temperature remains at a set point.  We follow this analysis with steepest-descent calculations to determine numerically the optimal heat source to maintain the outlet temperature at a fixed point.  We do this for both the case of small-amplitude inlet perturbations (linear theory, Section <ref>), as well as large-amplitude perturbations, which involves a non-linear theory (Section <ref>).    Discussion and concluding remarks are presented in Section <ref>.





 §.§ Scope of the work


The three-equation model has been formulated with disk driers in mind.  The simplified mathematical model is more generic.  Its main purpose is a test-bed for the optimal-control theory and a means of performing rigorous validation tests on our numerical methods.  Sections <ref>–<ref> on the simplified mathematical model also help us to formulate the optimal-control theory in a clear and transparent manner.  Then, it is possible to describe the optimal control theory for the three-equation model by simple analogy.  The main variable targeted for control in the present work is the outlet temperature.  In practical settings, it is important to maintain this variable near a set-point for food-hygiene and quality-control reasons.  Clearly, in an industrial drier, the outlet moisture is a key variable whose value needs to be controlled.  In the present work, this variable is maintained at the equilibrium value by the choice of design parameters – in particular, the drier length and the drier residence time.  Introducing a second manipulated variable (in addition to the heat source) with a view to gaining more control over the outlet moisture level will be the subject for future work, for which the mathematical framework introduced herein can be used.







§ SIMPLIFIED MATHEMATICAL MODEL



We look at a simplified mathematical model of drying.  This model involves only the product temperature, and reads as follows:




    ∂ T/∂ t+u_0∂ T/∂ x=k[q(t)-T],     x∈ (0,].

Here, u_0>0 is a constant velocity (dimensions: [Length]/[Time]), and k>0 is a rate (dimensions: 1/[Time]).  The quantity q(t) (with dimensions of temperature) represents the temperature of the surroundings, this may be manipulated to bring T(x,t) to a desired value at the outlet x= (i.e. (t)=T(,t)).


Hence, the  temperature at the inlet varies in time according to:

    T(0,t)=(t),


where (t) is a time-dependent inlet condition, with t≥ 0.  Furthermore, at t=0, the temperature is determined by an initial condition:


    T(x,0)=(x),

with x≥ 0.  The solution is continuous if (0)=(0).  





 §.§ Analytical and Numerical Solutions


Equation (<ref>) admits an analytical solution, which may be validated by direct computation:



    T(x,t)=(x-u_0 t)^-kt+(t),   u_0 t<x,
    [(t-x/u_0)-(t-x/u_0)]^-kx/u_0+(t),   u_0 t>x,

where 



    (t)=^-kt∫_0^t ^kt'[k q(t')]  t'.



Equation (<ref>) clearly highlights the hyperbolic nature of the  model in Equation (<ref>): at early times (u_0 t<x), information from the initial condition dominates, and is carried downstream by the constant velocity u_0.  At later times (u_0 t>x), information from the inlet dominates, and is carried downstream by the flow.  




 Equation (<ref>) can also be described numerically.  For this purpose,  the equation is discretized in space and time using a second-order upwind scheme:


    T_i^n+1-T_i/Δ t+u_0(3T_i^n-4T_i-1^n+T_i-2^n/2Δ x)=
    k[q(t^n)-T_i^n].

This is an explicit scheme, as such it is constrained by a CFL condition, u_0Δ t<Δ x.  The numerical approach is useful in caess where analytical solutions are not available, for instance, in the three-equation model.







 §.§ Optimal Control – Special Analytical Solution


For the Simple Mathematical Model, an analytical solution for the optimal control q(t) can be constructed, which maintains the outlet temperature T(,t) at a set-point T_*.  To see this,  we simply set T(,t)=T_*
in the analytical solution in Equation (<ref>).  This gives:



    T_*=(-u_0 t)^-kt+(t),     u_0 t<,

and


    T_*=[(t-/u_0)-(t-/u_0)]^-k/u_0+(t),     u_0 t > .



Multiply Equation (<ref>) by ^kt and differentiate with respect to time to obtain:



    T_*+u_0/k'(-u_0 t)^-kt+q(0)^-kt=q(t),     t</u_0.

A similar treatment of Equation (<ref>) gives:


    T_*-(t-/u_0)^-k /u_0-1/k'(t-/u_0)^-k/u_0
    
    +q(t-/u_0)^-k/u_0=q(t),     t>/u_0.




Equations (<ref>) contain an unknown quantity q(0).  However, in the case where the initial temperature, the mean temperature T_0, and the set-point temperature T_* are all the same (=T_0=T_*), q(0) can be taken as zero.  Then, Equation (<ref>) has a closed-form solution, and Equation (<ref>) can be solved recursively.  Notice that the optimal control has a possible jump discontinuity at the characteristic time t_0=/u_0, with


    q(t_0^-)   =   T_*+u_0/k'(0)^-kt_0,
    
    q(t_0^+)   =   T_*-(0)^-kt_0-1/k'(0)^-kt_0.





	
		
		
	






	
		
		
	








 §.§ Optimal Control – limit of large k/u_0


Assume a bounded control, |q(t)|≤ M, for all t∈[0,∞), and a bounded inlet condition |(t)|≤ M' for all t∈[0,∞).  Then, (t) is bounded by (t)≤ kM(1-^-kt).  Hence, in the limit of large k/u_0, the analytical solution (<ref>) simplifies:


    T(,t)=(-u_0 t)^-kt+(t),    t</u_0,
    (t),    t>/u_0,

Hence, the influence of the inlet at late times t>/u_0 is `washed away' by the exponential factor ^-k/u_0.  Thus, at late times, the outlet temperature can be manipulated by choice of (t).  Simply taking q(t)=T_* gives (t)=T_*(1-^-kt) and hence,


    T(,t)=T_*(1-^-kt),     t>/u_0,

Thus, the outlet condition at late times is controlled simply by taking q(t)=T_*, provided k/u_0≫ 1.





§ GENERAL OPTIMAL CONTROL THEORY FOR THE SIMPLIFIED MATHEMATICAL MODEL




To extend our method of control to cases where analytical solutions are not possible, we formulate a more general theory of optimal control based on a steepest-descent technique.

Hence, we introduce the penalty function


    J(T,q)=12∫_0^τ[T(,t)-T_*]^2 t.


Here, τ is the final simulation time.
The aim of this section is to develop a general framework for computing the optimal control q(t) at each point in time  t∈ [0,τ], such that J is minimized. 


For these purposes, we seek to minimize the functional J subject to the constraint that T(x,t) satisfies the mathematical model (<ref>).  Hence,   we introduce a
 a constrained functional:


    (T,q,ψ)=12∫_0^τ[T(,t)-T_*]^2 t+
    ∫_0^τ t∫_0^ x
    [∂ T/∂ t+u_0∂ T/∂ x
    -k(q(t)-T)]ψ.

Here, ψ is a Lagrange multiplier.  The physical interpretation fo ψ is that it enforces the constraint that T(x,t)  satisfies the model equations (in this context, Equation (<ref>)).   We follow standard practice in developing an equation for ψ.  However, we note that in some applications, the vanishing of the Lagrange multiplier engenders the so-called Lagrange catastrophe, which may cause the constraints to be violated <cit.>.  In this case, the semi-inverse method can be used to develop the appropriate functional to which a variational principle can be applied <cit.>.

We next compute δ, which is the variation in  with respect to variations in T, and variations in q (denoted by δ q).  Note that δ q depends only on time, not on space.  We compute δ in a straightforward fashion and obtain:



    δ=∫_0^τ[T(,t)-T_*]δ T(,t) t
    -∫_0^τ t∫_0^ kψδ q  x
    
    +∫_0^τ∫_0^[-∂ψ/∂ t-u_0∂ψ/∂ x+kψ]δ T  x+
    ∫_0^[ψδ T]_t=τ x+∫_0^τ[u_0ψδ T]_x= t.

Here, in computing the boundary terms (which arise from using integration by parts), we have assumed that δ T=0 at the inlet and also, δ T=0 at t=0.  This makes sense, since inlet/initial conditions are fixed, so there is no variation in the solution at these points.  Based on this expression for δ, we identify the adjoint problem:





    -∂ψ/∂ t-u_0∂ψ/∂ x+kψ=0,     x∈ [0,).

We impose the outlet condition 


    ψ(,t)=-(1/u_0)[T(,t)-T_*],

and we further impose the terminal condition 


    ψ(x,τ)=0.



With the adjoint problem defined in this way, the variation in δ simplifies dramatically:



    δ=-∫_0^T  t  δ q∫_0^ kψ(x,t) x.

These calculations now provide a numerical method for the computation of the optimal control q(t):











Operations 2–4 correspond to one iteration in a cycle.  Label the iterations by n, with n=0,1,2,⋯.  By carrying out repeated iterations, we obtain:


    δ=^n+1-^n=-λ^n∫_0^τ[∫_0^(kψ^n)  x]^2 t,

for sufficiently small λ^n, and correspondingly,


    q^n+1(t)=q^n(t)+λ^n ∫_0^ℓ (kψ^n) x.

Hence, δ<0 and the iterative process will converge, such that δ→ 0 as n→∞, such that  
becomes stationary, and the optimal control for q(t) is obtained.  Indeed, this algorithm is essentially a  steepest-descent method to compute the minimum of the functional , where the search direction is given by:


    -δ/δ q(t)=∫_0^ kψ(x,t) x.


With the steepest-descent method identified in this manner, it is possible to use standard optimization techniques to compute the optimal value of λ^n at each iteration cycle, this is addressed later in the section.


Finally, we notice that Equation (<ref>) is reminiscent of the Variational Iteration Method for solving non-linear ordinary differential equations <cit.>.  The Variational Iteration Method has also been used to solve problems in optimal control theory for ordinary differential equations <cit.> and also, partial differential equations <cit.>.  Hence, the Variational Iteration Method may be applicable to the present optimal-control problem as dwell, and would be complementary to the numerical method developed herein.





 §.§ Mathematical Setting


We place the foregoing calculations in a more mathematical context as follows.  We let T(x,t)∈ Y and q(t)∈ U.  We let e(T,q)=∂_t T+u∂_x T-k(q-T)∈ Z.  As e is linear in T and q, this can also be identified with a linear problem,


    e(T,q)=AT+Bq,     A∈ℒ(Y,Z).


Here, Y and Z are appropriate Banach spaces, ℒ(Y,Z) is the space of linear operators from X to Y, and U is an appropriate Hilbert space.  Hence, the weak form of the simplified mathematical model reads:


    e(T,q)=0 in Z,

where e: Y× U→ Z is defined with respect to its action on test functions:


    ⟨ψ,e(T,q)⟩_Z^*,Z=
    ∫_0^τ x∫_0^ (-T∂_t ψ)  x
    
    +∫_0^τ x∫_0^ u(-T∂_x ψ)  x
    -∫_0^τ x∫_0^ k(T-q)ψ  x
    
    +∫_0^τ u[T(,t)ψ(,t)-T(0,t)ψ(0,t)] t
    +∫_0^[T(x,τ)ψ(x,τ)-T(x,0)ψ(x,0)] x,

for all test functions ψ∈ Z^*.

The Lagrangian of the system is L:Y× U× Z^*→ℝ, where


    L(T,q,ψ)=J(T,u)+⟨ψ,e(T,q)⟩_Z^* Z

which we recognize as Equation (<ref>).  The first-order optimality conditions in the direction (δ T,δψ,δ u) are:



    ⟨ J_T(T,q),δ T⟩_Y^*,Y+⟨ψ,e_T(T,q)δ T⟩_Z^*,Z   =   0,
    ⟨δψ,e_T(T,q)δ T⟩_Z^*,Z   =   0,
    ⟨ J_q(T,q),δ q⟩_U^*,U+⟨ψ,e_q(T,q)δ q⟩_Z^*,Z   =   0;


these conditions must hold for all (admissible) directions δ T∈ Y, δψ∈ Z, and δ q∈ U.  We recognize Equation (<ref>) as the adjoint problem (<ref>), Equation (<ref>) as the weak formulation of the simplified mathematical model (Equation (<ref>)), and Equation (<ref>) as the optimality condition δ J/δ q=0 (cf. Equation (<ref>)).

Furthermore, we look at solutions to e(T,q)=0, such that T depends implicitly on the forcing term, T=T(q).  Hence, we look at the reduced cost function,


    J(q)=J(T(q),q)=J(T(q),q)+⟨ψ,e(T(q),q)⟩_Z^*,Z=L(T(q),q,ψ).


Differentiating this expression in a direction δ q∈ U, we obtain:


    ⟨J'(q),δ q⟩_U^*,U=⟨ L_T(T(q),q,ψ),∂ T/∂ qδ q⟩_Y^* Y+⟨ L_q(T(q),q,ψ),δ q⟩_U^*,U.

We now choose ψ=ψ(q) such that:


    L_T(T(q),q,ψ)=0

(this is nothing but the adjoint equation (<ref>)).  
In fact,


    ⟨ L_T(T,q,ψ),d⟩_Y^*,Y   =   ⟨ J_T(T,q),d⟩_Y^*,Y+⟨ψ,e_T(T,q)d⟩_Z^*,Z,
       =   ⟨ J_T(T,q)+e_T(T,q)^*ψ,d⟩_Z^*,Z.

Furthermore, we now have:


    ⟨J'(q),δ q⟩_U^*,U   =   ⟨ L_q(T(q),q,ψ(q)),δ q⟩_U^*,U,
       =   ⟨ J_q(T(q),q),δ q⟩_U^*,U+⟨ψ(q),e_q(T(q),q),δ q⟩_U^*,U,
       =   ⟨ J_q(T(q),q),δ q⟩_U^*,U+⟨ e_q(T(q),q)^*ψ(q),δ q⟩_U^*,U,

which is true for all admissible perturbations δ q∈ U.  As U is a Hilbert space, the Riesz representation theorem applies, and we can identify the vector J'(q)∈ U^*=U:


    J'(q)=J_q(T(q),q)+e_q(T(q),q)^*ψ(q),

which is precisely Equation (<ref>).

In Section <ref>, we have already established (by direct calculation) that the  optimization problem


    min_(T,q)∈ Y× UJ(T,q),    subject to  AT+Bq=0,     q∈ U,     T∈ Y,

has a minimum – specifically, that min_(T,q)∈ Y× UJ(T,q)=0.   We also look here at the more general case,


    min_(T,q)∈ Y× UJ(T,q)+12αq_U^2,subject to  AT+Bq=0,
     u∈ U_ad⊂ U,     T∈ Y_ad⊂ Y,

Then, the following result holds<cit.>:

  (Hinze et al.)  Suppose the following statements are true:



  * α≥ 0, U_ad⊂ U is convex, and in the case of α=0, bounded;

  * Y_ad⊂ Y is convex and closed, such that AT+Bq=0 for at least one set of functions T∈ Y_ad and q∈ U_ad;

  * A∈ℒ(Y,Z) has a bounded inverse.

Then, the problem (<ref>) has an optimal solution (T_opt,q_opt).  If α>0, then the solution is unique.


For the present case, the form of J suggests the choice U=L^2([0,τ]).  
We limit ourselves to bounded controls, such that q∈B(0,r)=U_ad⊂ U.  We also take Y_ad=Y, 
and A=∂_t+u_0∂_x+k.   The explicit solution of AT=kq given in Equation (<ref>) establishes the that A has a bounded inverse.  Thus, the conditions of the Theorem are satisfied, and a unique optimal solution exists. 

Beyond the Simple Mathematical Model, we look at a detailed three-equation model in Sections <ref>–<ref>.  We look at the linearized three-equation model in  Section <ref>, to which Theorem <ref> again applies.  In this case, Fourier / Laplace theory can be used to construct an explicit solution for T(q), which establishes the existence of a bounded inverse operator for the problem.  We furthermore look at the nonlinear version of the model in Section <ref>.  This takes the form of a set of coupled semi-linear parabolic equations; the equations can be reduced to a single semi-linear equation for temperature T, with a continuous non-linear term.  In this case, Reference <cit.> establishes the existence of an optimal control.  Other works which look at nonlinear differential equations from the point of view of variational calculus include References <cit.>, and <cit.>.







 §.§ Barzilai–Borwein Method




In the present work, we are interested in minimizing a cost function min_x∈ X J(x).  In practice, this amounts to solving an equation of the type


    G(x)=0,


in this context, G(x) is the gradient of the cost function J.  Here, X is an appropriate Hilbert space, with inner product ⟨·,·⟩_X, and with norm x_2^2=⟨ x,x⟩_X.  The generalized Newton method is used to find the solution x of Equation (<ref>)<cit.>:




There are a number of methods for choosing the operator M_k.  In the steepest-descent method, M_k is chosen to be a diagonal matrix, M_k=λ_k^-1𝕀, thus, the sequence of steps {x_0,x_1,⋯,x_k,⋯} follows the path of steepest descent towards the minimum of J(x).   The length of the descent step λ_k is chosen such that F(x_k+1)<F(x_k), and the method is known to possess linear convergence.



The classical Newton method involves choosing M_k to be the derivative of G.  In this case, convergence of the iterative method is super-linear, provided the derivative of G is a continuously invertible linear operator on X<cit.>.  However, often it is expensive to compute the derivative of G, in which case the derivative of G may be approximated by the Secant Method:


    d_k=x_k-x_k-1,     g_k=G(x_k)-G(x_k-1),     M_k d_k=g_k,    (Secant Method).


The Barzilai–Borwein method combines the ideas of the Steepest-Descent method and the Secant Method.
Hence, the matrix M_k in the Secant Method is approximated as a diagonal matrix:



    M_k=1/α_k𝕀.



The equation M_k d_k=g_k  is then solved in the least-squares sense:


    α_k=arg min_α>01/αd_k-g_k_2^2,

which gives


    α_k=d_k_2^2/⟨ d_k,g_k⟩_X.

Then, the update step (<ref>) in the Barzilai–Borwein method becomes:


    x_k+1=x_k-α_k g(x_k).




Although the Barzilai–Borwein method does not decrease the objective function monotonically,  extensive numerical experiments show that it performs
much better than a standard Steepest-Descent method<cit.>.  In the present work, we use the Barzilai–Borwein method exclusively, and monitor the convergence of the method (and the corresponding decrease in the objective function) carefully.








 §.§ Sample Numerical Results



We use the Barzilai–Borwein method to compute the optimal control q(t) for the model problem (<ref>), with =T_*=100^∘ C, time in minutes, and the other parameters given in Table <ref>.  The parameter values in Table <ref> are not from a real industrial drier but rather, are representative of the timescales observed in such a drier.  A detailed set of representative parameters for an industrial drier is given in Section <ref>, in Table <ref>.


The inlet forcing is given by:


    T_inlet(t)=[100+10sin(2π t)]^∘C.








 A comparison between the numerical optimal control and the analytical  one  is shown in Figure <ref>. 



Near the  characteristic time t_0=/u_0, there is some discrepancy between the two approaches.  This comes from the underlying discretization of the model partial differential equation (i.e. Equation (<ref>)), which through numerical diffusion, introduces some extra smoothness into the numerical solution.  Otherwise, there is excellent agreement between the two methods, and this validates the numerical approach to computing the optimal control.




Convergence of the Barzilai–Borwein method is shown in Figure <ref>.  The convergence is non-monotone, however, each `spike' in the value of the cost function J is followed by a subsequent reduction in J, which leads to an overall reduction in J to close to a value  of 10^-10 after 1000 iterations.







 §.§ Inlet Forcing


In the foregoing results, we have used sinusoidal inlet forcing (Equation (<ref>)).  Often, in the food-processing industry, the inlet temperature is characterized by a stochastic process, and as such the inlet temperature consists of a mean value, plus a random perturbation.  The random perturbations can be quite complex (i.e. more complicated than a simple Brownian process), and exhibit correlations over a timescale of hours.


  We have therefore looked at using time series of temperature obtained from sensor data from real industrial driers as an inlet forcing term, both for the simplified mathematical model and the three-equation model.  The two types of inlet forcing  produce optimal controls with the same type of behaviour.  Hence, for the purpose of introducing and validating the present mathematical framework, results are presented for sinusoidal forcing throughout this article.  The applicability of the results to stochastic inlet forcing is discussed in Section <ref>.







§ INDUSTRIAL DRIER: THEORETICAL FORMULATION



The operating principle of an industrial disk drier is shown schematically in  Figure <ref>.  Product is fed in to the drier at the inlet, the product consists of solids and water in the liquid phase, at a given inlet temperature .  The product is then carried through the drier at constant velocity u_0 and subjected to a heat source (direct contact, via heated disks in the case of the disk drier).   The product is assumed to be well mixed, such that the dry solid component and the water component are at a common temperature T.  The product is exposed to a moist gas phase, the heat source heats up the product and causes the water in the product to evaporate into the moist gas.  A continuous stream of air enters the drier at some location; also, the moist gas is vented from the drier.  There is therefore a continuous flow of moist gas through the drier, the flow rate is assumed to be large, such that the moist gas phase is turbulent and well mixed  and hence, the moist gas phase is treated as homogeneous and steady.


The dry solid volume fraction is denoted by χ_s, the volume fraction of the water in the solid is denoted by χ_l, and the volume fraction of the moist gas is χ_g, with the consistency condition χ_s+χ_l+χ_g=1.


 

Assuming that the conditions of the drier are homogeneous in the plane perpendicular to the page, only the dimension along the length of the drier (x-direction) is important.  Hence, we identify the  three variables as determining the state of the drier at any instant: ϵ_s=ρ_s χ_s, ϵ_l=ρ_l χ_l, and T, the product temperature.  The internal energy e per unit mass of product is a function of these variables, e=e(ϵ_s,ϵ_l,T).  Conservation of mass then yields:



    ∂ϵ_s/∂ t+∂/∂ x(u_0 ϵ_s)   =   0,
    ∂ϵ_w/∂ t+∂/∂ x(u_0 ϵ_w)   =   -ṁ.





Here, ṁ is the drying rate, with dimensions of kg/(m^3·s), and is a general function of ϵ_s,ϵ_l, and T.  

The density of the product phase is denoted by ρ=ϵ_s+ϵ_l, this satisfies the equation


    ∂ρ/∂ t+∂/∂ x(u_0 ρ)=-ṁ.


The energy balance equation is formulated in terms of e, the internal energy per unit mass in the product stream.  This is given by:


    e=[∑_i=s,l c_p,i(ϵ_i/ρ)](T-).

Here, T is the temperature of the product,  is a reference temperature, and the ratios ϵ_i/ρ play the role of mass fractions.  Also, c_p,i is the specific heat at constant pressure (per unit mass basis) in the i^th component of the product mixture. Thus, the equation for energy balance in the product stream reads:


    ∂ρ e/∂ t+∂/∂ x(u_0 ρe)=q̇-ṁ h_l.


Here, q̇ is the heat transfer into the product (units: J/(m^3·s)), and h_l is the latent  heat of evaporation, per unit mass.   Equation (<ref>) can be further simplified to read:


    ∂ T/∂ t+u_0∂ T/∂ x=q̇-ṁ[h_l-c_p,w(T-)]/∑_i=s,lc_p,iϵ_i,


which is an advection equation for temperature.
Thus, there are three state variables: the densities ϵ_w and ϵ_l, and the product temperature T.

Two laws  / closure models are required (one for q̇ and one for ṁ).   The three  model equations (<ref>), (<ref>), and (<ref>) transfer inputs at the inlet to outputs at the outlet.  











































 



































 §.§ Drying rate


The functional form of the drying rate has two extreme cases.  In the initial drying stage, the drying rate is limited by the ability of the moist gas to take up evaporated water from the product.  This is the so-called constant drying-rate regime.  In this regime, the drying rate reads<cit.>:


    ṁ=k_c a_V(Y_*(T)-Y_air).

 Here, k_c is a constant with units of kg/(m^2·s), a_V is the surface area of the product in contact with the moist gas, per unit volume, Y_* is the absolute humidity of the air at saturation, at the product temperature T.  Also, Y_air is the absolute humidity of the air.  Here, the absolute humidities are expressed on a dry-gas basis, that is, Y=M_wv/M_da, where M_wv is the mass of water vapour in a parcel of moist gas, and M_da is the amount of dry gas in that same parcel.  













Beyond a certain critical moisture level in the solid, there is no longer sufficient water at the product surface to enable direct evaporation.  The critical moisture level is given in terms of the  absolute moisture content of the product on a dry-solid basis, X=ϵ_w/ϵ_s, the critical level is denoted by X_c.

 Below this level, the drying is diffusion-dominated, as it relies on diffusive processes to bring water from the interior of the product to the surface, where evaporation can then occur.  This is the falling drying-rate regime, in this regime the drying rate is given by the Lewis equation<cit.>,


    ṁ=k_f ϵ_s(X-X_*).

Here, k_f is a rate coefficient (units: s^-1), and X_* is the equilibrium moisture content of the solid.  In some applications, k_f is known to be a function of the product temperature.  Equation (<ref>) can be rewritten as:


    ṁ=k_f (ϵ_l-X_*ϵ_w);

Summarizing, the expression for the drying rate which accounts for the two different regimes is:


    ṁ=k_c a_V(Y_*(T)-Y_air),   X>X_c,
    
                        k_f (ϵ_l-X_*ϵ_w),   X≤ X_c.

where k_c and k_f are the two rate coefficients which are not known a priori.

In the present article, we are concerned mostly with the optimization of drying, which requires functional differentiation of the drying rate.  Hence, the discontinuous drying rate in Equation (<ref>) is not suitable.  In practice, the drying-rate curve (which encodes all information about ṁ) exhibits a smooth transition from the constant-rate regime to the falling-rate regime.  Therefore, a more appropriate functional form of the drying rate would be to have a smooth interpolation between the two regimes in Equation (<ref>).  However, for the purposes of showcasing the optimization techniques which we develop in this article, we instead focus  exclusively on the falling drying-rate regime, and use the expression ṁ=k_f(ϵ_l-X_*ϵ_w), with the assumption that both k_f and X_* are constants.



 §.§ Heating


We consider a constant amount of thermal energy delivered to the product along the length of the drier, this can be achieved for instance by indirect contact of the product with a closed-loop flow of steam.  The thermal energy per unit time, per unit volume is denoted by q̇.  Thus, the total power is


    P=A_×∫_0^q̇  x,

where A_× is the cross-sectional area of the drier.  It is assumed that the power is delivered to the product evenly along the length of the drier, thus, we can write 


    q̇=P/A_×.













§ INDUSTRIAL DRIER: EQUILIBRIUM SOLUTIONS, LINEARIZATION



In this section we consider again the three-equation model of the industrial drier, which we summarize here as:




    ∂ϵ_s/∂ t+∂/∂ x(u_0 ϵ_s)   =   0,
    ∂ϵ_w/∂ t+∂/∂ x(u_0 ϵ_w)   =   -ṁ,
    ∂ T/∂ t+u_0∂ T/∂ x   =   q̇-ṁ[h_l-c_p,w(T-)]/∑_i=s,lc_p,iϵ_i_=H,



where ṁ=k_f(ϵ_l-X_*ϵ_s).
We focus first of all on developing an equilibrium solution of this model, and then develop a linear stability analysis, outlining the conditions under which the equilibrium solution is stable to small perturbations at the inlet.  







 §.§ Equilibrium Solution



We consider Equation (<ref>) with fixed inlet conditions ϵ_s(0,t)=ϵ_s0, ϵ_w(0,t)=ϵ_w0, and T(0,t)=T_0, where ϵ_s0, ϵ_w0, and T_0 are constants.  Under these inlet conditions, we seek an equilibrium solution of Equation (<ref>). We immediately obtain ϵ_s(x)=ϵ_s0, the constant solution.  The other two equations in the set simplify:


    u_0ϵ_w/  x   =   -k_f(ϵ_w-X_*ϵ_s),
    
    u_0  T/  x   =   H,



We continue with the assumption that q̇ is constant, and equal to q̇=P/(A_×), where P is the power.








An explicit expression for ϵ_w(x) can be obtained using the integrating-factor technique:



    ϵ_w(x)=ϵ_w0^-k_f x/u_0+ϵ_s0X_*(1-^-k_f x/u_0).


Then, back-substitution into Equation (<ref>) reduces the three-equation model down to a single equation for temperature.
However, for our purposes, it is more straightforward to solve the three-equation model numerically.   Hence, numerical solutions to the three-equation model are shown in Figure <ref>.








Here, we have used the constant value k_f=1/(5 min), based on a typical drying time for solids in a dryer.  We further use the typical parameter values in 
Table <ref> 

 (the velocity in Table 2 corresponds to a residence time of 30 mins, T_R=/u_0).





 §.§ Péclet Number


Based on Reference <cit.>, we take a value of k=0.6 W/(kg·K), the same for each component in the product (liquid and solid).  The thermal diffusivity is therefore k/∑_i=s,lc_p,iϵ_i, and the Péclet number is:


    Pe=u_0 ℓ/k/∑_i=s,lc_p,iϵ_i.

Based on the values in 

 Table <ref>, 
and using the values at the inlet for ϵ_i, the Péclet number is 8172.  The large value of the Péclet number justifies the neglect of diffusion on the three-equation model.  Furthermore, the absence of sharp temperature gradients in the numerical solutions suggests that there are no `boundary-layer effects', i.e. small regions where the diffusion term is important.  These results justify the choice of a hyperbolic model in the present work.



 §.§ Linear Stability Analysis


We consider a time-dependent inlet condition


    ϵ_s(0,t)=ϵ_s,0+δϵ_s(t),    ϵ_l(0,t)=ϵ_l,0+δϵ_l(t),    
    T(0,t)=T_0+δ T(t),

where δϵ_s(t), δϵ_l(t), and δ T are small-amplitude time-dependent fluctuations.  Hence, the solution of the three-equation model (<ref>) can be decomposed as:


    ϵ_s(x,t)=ϵ_s,0+δϵ_s(x,t),    ϵ_l(x,t)=ϵ_l,eq(x)+δϵ_s(x,t),    
    T(x,t)=T_eq(x)+δ T(x,t),

where ϵ_l,eq(x) and T_eq(x) denote the equilibrium solution.  We substitute Equation (<ref>) into the three-equation model (<ref>) and linearize around the equilibrium solution to obtain:


    ∂/∂ t([ δϵ_s; δϵ_l;  δ T ])
    +u_0∂/∂ x([ δϵ_s; δϵ_l;  δ T ])
    =𝕁([ δϵ_s; δϵ_l;  δ T ]),

where 𝕁 is the Jacobian matrix


    𝕁=([        0        0        0; -∂ṁ/∂ϵ_s -∂ṁ/∂ϵ_l  -∂ṁ/∂ T; ∂ H/∂ϵ_s ∂ H/∂ϵ_l  ∂ H/∂ T;          ])


We decompose the inlet perturbations into the harmonic components, using


    δ T(t)=1/2π∫_-∞^∞δ T_ω^ω tω,

and similarly for the inlet densities.  The perturbations in the body of the drier (i.e. x∈(0,]) can also be assumed to have the same decomposition.  Thus, the linear problem (<ref>) can be re-written in Fourier space as:


    u_0∂/∂ x([ δϵ_s,ω; δϵ_l,ω;  δ T_ω ])
    =(-ω𝕀+𝕁)([ δϵ_s,ω; δϵ_l,ω;  δ T_ω ])


The solution of Equation (<ref>) can be formally written as:


    ([ δϵ_s,ω(x); δϵ_l,ω(x);  δ T_ω(x) ])=
    ^[-ω x + ∫_0^x 𝕁 x]/u_0([ δϵ_s,ω(0); δϵ_l,ω(0);  δ T_ω(0) ]).


Hence, the question as to whether initial small-amplitude disturbances at the inlet decay to zero as they move downstream is determined by the eigenvalues of ∫_0^x𝕁 x.

As this matrix has a row of zeros, one of the eigenvalues is itself zero, therefore the characteristic equation simplifies:


    |[ -∫_0^x ∂ṁ/∂ϵ_l x-λ    -∫_0^x ∂ṁ/∂ T x;    ∫_0^x∂ H/∂ϵ_l x   ∫_0^x∂ H/∂ T x-λ ]|=0.


We first look at the case where k_f and X_* are independent of T, in this case there is a single negative eigenvalue and a single positive eigenvalue; the positive eigenvalue is given by:







    λ=∫_0^x∂ H/∂ T x=∫_0^x (ṁ c_p,w/c_pwϵ_w+c_p,lϵ_l) x≥ 0.


Hence,  small-amplitude inlet disturbances are amplified as they are convected downstream to the outlet.  





However, as the drier has finite length, the degree to which the inlet disturbances are amplified is also limited by the growth factor ^max_[0,]λ/u_0






This result applies to all inlet forcing frequencies ω.  Thus, the instability is reminiscent of a convective instability<cit.> – effectively, an amplifier, which amplifies all inlet disturbances at the same rate as they are convected downstream.  This is in contrast to an absolute instability, which is present in systems with a natural frequency.  It also contrasts to the simplified mathematical model in Section <ref>, which is a stable system, in the sense that disturbances at the inlet are damped to zero as they propagate to the outlet.



As the perturbations away from the inlet are bounded, it is possible to introduce a proportional (bounded) perturbation in the control variable q̇ to bring the outlet temperature T(,t) to a desired value.  To show this, we consider again the linear problem (<ref>), however, we now consider further the possibility that the heat source q̇ is also subject to a perturbation, which the operator can control.  Hence, we consider q̇=q̇_0+δq̇, where q̇_0 is a constant and equal to (P/A_×), and δq̇ depends only on time.




    ∂/∂ t([ δϵ_s; δϵ_l;  δ T ])
    +u_0∂/∂ x([ δϵ_s; δϵ_l;  δ T ])
    =𝕁([ δϵ_s; δϵ_l;  δ T ])+δq̇(t) ([    0;    0; η(x) ]),

where η(x)=[c_p,lϵ_l(x)+c_p,sϵ_s,0]^-1.  In Fourier space, this reads:


    u_0∂/∂ x([ δϵ_s,ω; δϵ_l,ω;  δ T_ω ])
    =(-ω𝕀+𝕁)([ δϵ_s,ω; δϵ_l,ω;  δ T_ω ])
    +δq̇_ω([    0;    0; η(x) ]).






















                         
												 
												
												
												
















												
												
												
												
												




												
												
												
												
												




												
												
												
												
												







We look again at Equation (<ref>), with the same restrictions as before that k_f and X_* are constants.  Then, there is an explicit solution for δϵ_s,ω(x) and δϵ_l,ω(x):


    δϵ_s,ω(x)   =   δϵ_s,ω(0)^-ω x/u_0,
    δϵ_l,ω(x)   =   δϵ_l,ω(0)^-ω x/u_0-k_fx/u_0+X_*δϵ_s,ω(0)^-ω x/u_0(1-^-k_fx/u_0).

Thus, inlet oscillations are propagated undamped to the outlet, this is due to 

 the zero eigenvalue of the Jacobian 𝕁.

Correspondingly, the equation for the temperature fluctuations reads:


    / xδ T_ω+(ω - ∂ H/∂ T/u_0)δ T_ω=
    1/u_0(∂ H/∂ϵ_sδϵ_s,ω+∂ H/∂ϵ_lδϵ_s,ω)_=ρ(x)+δq̇_ωη(x)/u_0,

with solution


    δ T_ω(x)=δ T_ω(0)^-ω x/u_0+λ_+(x)/u_0
    
    +^-ω x/u_0+λ_+(x)/u_0∫_0^x [ρ(x')/u_0]^ω x'/u_0-λ_+(x')/u_0 x'
    
    +δq̇_ω^-ω x/u_0+λ_+(x)/u_0∫_0^ [η(x')/u_0]^ω x'/u_0-λ_+(x')/u_0 x'.

Here, λ_+(x) is the positive eigenvalue in Equation (<ref>).
Hence, in order to suppress the temperature fluctuations at the outlet, such that δ T_ω()=0, it suffices to apply the control:


    δq̇_ω=-δ T_ω(0)-∫_0^ [ρ(x')/u_0]^ω x'/u_0-λ_+(x')/u_0 x'/∫_0^ [η(x')/u_0]^ω x'/u_0-λ_+(x')/u_0 x',

which exists, provided the denominator is non-zero.  

To account for the hyperbolic nature of the system, great care must be taken in reconstructing the solution of δ T(x,t) in the spatio-temporal domain, from the data given in δ T_ω(x).  In particular, the reconstruction needs to take account of the jump which occurs at x=u_0 t, where the solution δ T(x,t) transitions from being dominated by information at the initial time, to information from the inlet.   This was straightforward in the case of the simplified model (Section <ref>), but is considerably more difficult here because the forcing term δq̇(t)η(x) depends on space and time.  We do not pursue this approach any further here, and instead pass over to a numerical approach.





 §.§ Numerical Results








                          







	
		
		
	








	
	  
		
		
	




We apply optimal control theory to the linear system (<ref>), with a view to keeping the outlet temperature fluctuation at zero, δ T(,t)=0.  We use the numerical technique developed in Section <ref>, for which the corresponding adjoint problem reads:

    -∂/∂ t([ ψ_s; ψ_l; ψ_T ])
    -u_0∂/∂ x([ ψ_s; ψ_l; ψ_T ])
    -𝕁^T([ ψ_s; ψ_l; ψ_T ])=0,

with terminal condition ψ(x,τ)=0 at the end-time τ and outlet conditions


    ψ_s(,t)=0,ψ_l(,t)=0,    ψ_T(,t)=-(1/u_0)δ T(,t).

The n^th guess for the optimal control is labelled by [δq̇(t)]^n, and the updated guess at the next iteration is given by


    [δq̇(t)]^n+1=[δq̇(t)]^n+λ∫_0^Ψ_T(x,t) η(x,t) x,

where η(x,t)=[c_p,lϵ_l(x,t)+c_p,sϵ_s(x,t)]^-1, and
where λ is a small parameter that is chosen using the Barzilai–Borwein method.

We look at  monochromatic inlet forcing




    δ T(t)=(5^∘C)sin(ω_F t),

where ω_F=2π/T_F, and the forcing period T_F is taken to be 8.5 min.  Also, the inlet solid and liquid density fluctuations are set to zero.

Numerical results are shown in Figure <ref> for 1,000 iterations of the steepest-descent algorithm, with a final residual value 0.0378, this can be reduced further by applying further iterations of the steepest-descent algorithm.

We  use N=200 gridpoints and Δ t=0.1 s, for converged numerical results.   The optimization is carried out over a four-hour time horizon (i.e. τ=4 hr).





We focus our attention on the behaviour of δq̇, shown in Figure <ref>(d).  This consists of a single dominant oscillation with frequency ω_F, i.e. to match the inlet forcing.  There is a second, long-time trend present, in the form of a `beat', as evidenced by the long-time periodic modulation of the optimal control δq̇ in the figure.  To understand this in more detail, we have computed the power spectrum of the optimal control, shown in Figure <ref>.

There is a single sharp peak in the data at ω≈ω_F (specifically, ω=0.01222) and a secondary peak nearby, this time at ω=0.01047.  
The secondary peak is due to the `ramp-up' in the forcing term which occurs at x=u_0 t, where the solution goes from initial-condition-dominated to inlet-dominated.  This is confirmed by looking at other (non-optimal) forcing terms δq̇ – only those with the `ramp-up' exhibit a secondary peak.



As these two peaks are close in frequency, this gives rise to a `beating' phenomenon: the beat frequency is ω_B=0.01222-0.01047 (no division by two), which corresponds to a beat period of 1 hours, this manifests itself in the real-space signal as the long-period modulation.  








§ INDUSTRIAL DRIER: NONLINEAR OPTIMAL CONTROL THEORY



We now apply optimal control theory to the nonlinear system (<ref>), with a view to keeping the outlet temperature fluctuation at a set point T(,t)=T_*.  We use the numerical technique developed in Section <ref>.  Hence, we solve the nonlinear problem forward in time with a guess q̇(t) for the optimal forcing, specifically, we solve:



    ∂ϵ_s/∂ t+∂/∂ x(u_0 ϵ_s)   =   0,
    ∂ϵ_w/∂ t+∂/∂ x(u_0 ϵ_w)   =   -ṁ,
    ∂ T/∂ t+u_0∂ T/∂ x   =   q̇-ṁ[h_l-c_p,w(T-)]/∑_i=s,lc_p,iϵ_i_=H.



Next, we solve the adjoint problem:

    -∂/∂ t([ ψ_s; ψ_l; ψ_T ])
    -u_0∂/∂ x([ ψ_s; ψ_l; ψ_T ])
    -𝕁^T([ ψ_s; ψ_l; ψ_T ])=0,

with terminal condition ψ(x,τ)=0 at the end-time τ and outlet conditions


    ψ_s(,t)=0,ψ_l(,t)=0,    ψ_T(,t)=-(1/u_0)[ T(,t)-T_*].

Here, the Jacobian 𝕁 has exactly the same form as in the linear stability analysis (Section <ref>), Hence,  although the forward problem is nonlinear, the adjoint problem is linear.

However,  the Jacobian must now be evaluated on the densities and temperature from the forward-problem, which vary in space and time.  Consequently, knowledge of the forward variables at all points in space and time is required for the adjoint computation – this means that the simulation time is much longer than before in the case of the fully linear problem.   These are  generic features of optimal control problems for partial differential equations.     As before, we label the n^th guess for the optimal control  by [q̇(t)]^n, and the updated guess at the next iteration is given by


    [q̇(t)]^n+1=[q̇(t)]^n+λ∫_0^ψ_T(x,t) η(x) x,

where λ is a small parameter that is chosen (again) using the Barzilai–Borwein method.  This approach is iterated until the outlet temperature is sufficiently close the set point, for all times.




 §.§ Results


For the initial condition of the forward problem, we use the equilibrium solution from Section <ref>, ϵ_s,0 (a constant), ϵ_l,eq(x), and T_eq(x).  The inlet conditions are then selected as:



    ϵ_s(0,t)=ϵ_s,0,
    ϵ_l(0,t)=ϵ_l,eq(0)[1+δαsin(ω_F t)],
    
    T(0,t)=T_eq(0)[1+δαsin(ω_F t)],




where ω_F is the forcing frequency, which we again take to be ω_F=2π/(8.5 mins).   We take the magnitude of the initial inlet perturbation δα to be 0.05.  Also, in order to avoid a discontinuous control at t=0, we also take the set point T_* to be T_eq(). 

Results are shown in Figure <ref>, for a simulation  out to τ=2 h.




These are qualitatively very similar to the previous linear results.  The optimal control q̇(t) keeps the outlet temperature at the set point T_*.  A power-spectrum analysis of q̇(t) reveals a single peak at the forcing frequency, as well as neighboring peaks, thus producing a `beat' pattern in the optimal forcing (Figure <ref>) – qualitatively the same as in the linear case.  The large peak in the power spectrum in Figure <ref> at ω=0 corresponds to the constant component, which is treated separately in the linear stability analysis.






We have also explored large-amplitude initial disturbances using the same initial condition as before (Equation (<ref>)), but with δα=0.2.  In this case, the optimal control has q(t)≤ 0 for some t, which corresponds to cooling the product, as opposed to heating.  As this is technically infeasible, we look at imposing a constraint q(t)≥ 0 for all t∈[0,τ], this can be achieved for instance by letting q(t)=(1/2)[θ(t)]^2, and performing the steepest-descent calculations using δ/δθ=θ(t) ∫_0^ψ_T(x,t)η(x) x.


Results are shown in Figure <ref>, and are now qualitatively different from the previous linear and small-amplitude nonlinear cases.  As we are now constrained to have q̇≥ 0, the set point T(,t)=T_* is not reached, instead, the difference between T(,t) and T_* is minimized in the root-mean-square sense.







However, the optimal control in Figure <ref>(c) still  exhibits the previous `beat' pattern.  




	
		
		
	








 §.§ Outlet Moisture Level


From Figures <ref>–<ref> (in particular, panel (b) in both figures), it can be seen that the outlet moisture level X(,t) remains close to the equilibrium moisture content X_*.  This is a simple consequence of our choice of parameters for the drier – the velocity u_0, and the length of the dryer have both been chosen with reference to the value of k_f such the outlet equilibrium moisture level X_eq(,t) is close to X_*.  Furthermore, inspection of the linear stability analysis in Section <ref> shows that the liquid's density fluctuations are damped by the exponential factor ^-k_fx/u_0 where again, the value of  means that significant damping of fluctuations occurs along the length of the channel.  Thus, the effective control of the outlet moisture levels can be thought of as byproduct of good design.  

In contrast, the attainment of the control T(,t)→ T_* is a consequence of active control whereby q̇(t) is selected according to the steepest-descent method.








§ DISCUSSION AND CONCLUSIONS



Summarizing, we have developed a mathematical framework for spatio-temporal control in industrial drying, starting from first principles, namely convection-type models for heat and mass transfer.    In contrast to previous models, ours take account of the temperature variations of product which can occur inside the drier itself, these can be significant.

We have applied optimal control theory to the models, using a Steepest-Descent Method with the Barzilai–Borewein formula for stepsize, to determine the optimal external heating required to maintain the outlet temperature close to a set-point-temperature.  For the three-equation model, and in the case of small-amplitude inlet forcing, the optimal control maintains the outlet temperature at the set point; for large amplitudes this is not possible, and the optimal control minimizes the root-mean-square difference between the outlet temperature and the set-point temperature.

The main aim of the present work is to control the outlet temperature of the product, using a time-dependent heat source  as a single control variable. The outlet moisture level has been maintained at the fixed, equilibrium value by virtue of a judicious choice of design parameters – including the residence time and the drier length.  In cases where these parameters cannot be chosen in this way, a second manipulated variable should be considered – for instance, a valve at the outlet, allowing for the product to stay in the drier for longer.  In this case, the general mathematical framework introduced in the present work still applies.


The optimal control theory requires advance knowledge of the inlet conditions. 

In practical settings, this may not be available; however, knowledge of the state of the plant obtained from the model will be supplemented by sensor data for temperature at various points along the drier.  In this case, a Kalman Filter can be used to correct the model, in order to keep the model predictions `close' to observations.  The Kalman Filter may include process noise – which may be appropriate for modelling inlet conditions that consist of a mean value plus fluctuations.  This approach can be supplemented with model predictive control – a key  methodology from control theory used for estimating the optimal control based on the given mathematical model, the sensor data coming from the plant, and the Kalman Filter.
Our work provides a theoretical basis for understanding how model predictive control would work in case of an industrial drier with temporal and spatial variations.  First, the model predictive control relies on a robust discretized numerical model – such as the ones developed herein.  Second, the optimal controls presented herein 
 still form a `best case scenario' against which the performance of an online model predictive control can be assessed.  Hence, the models and methods developed herein form a key validation step which may be developed further in a more practical setting.





 §.§ Acknowledgements




The author acknowledges insights from discussions with       James Mc Carthy of Cervos Ltd., as well as from Luis Carnevale and Miguel Bustamante.  The author finally acknowledges the work done by Odhrán Dooley in helping to compile the parameter values in Table <ref>.




 §.§ Funding


This work has been produced as part of ongoing work within the ThermaSMART network.
The ThermaSMART network has received funding from the European Union's Horizon 2020
research and innovation programme under the Marie Sklodowska–Curie grant agreement No.
778104.








10

mujumdar2006handbook
Arun S Mujumdar.
Handbook of industrial drying.
CRC press, 2006.

hernandez2014modeling
Andres Hernandez, Robin De Keyser, José Manrique, Jose Oliden, and William
  Ipanaqué.
Modeling and nonlinear model predictive control of a rotary disc
  dryer for fishmeal production.
In 2014 European Control Conference (ECC), pages 1819–1824.
  IEEE, 2014.

luz2010dynamic
Gianini Regina Luz, Wagner André dos Santos Conceição,
  Luiz Mário de Matos Jorge, Paulo Roberto Paraíso, and Cid
  Marcos Gonçalves Andrade.
Dynamic modeling and control of soybean meal drying in a direct
  rotary dryer.
Food and bioproducts processing, 88(2-3):90–98, 2010.

iguaz2003mathematical
A Iguaz, A Esnoz, G Martınez, A López, and P Vırseda.
Mathematical modelling and simulation for the drying process of
  vegetable wholesale by-products in a rotary dryer.
Journal of food engineering, 59(2-3):151–160, 2003.

evans1983introduction
Lawrence C Evans.
An introduction to mathematical optimal control theory version 0.2.
Lecture notes available at http://math. berkeley. edu/
  evans/control. course. pdf, 1983.

hinze2008optimization
Michael Hinze, René Pinnau, Michael Ulbrich, and Stefan Ulbrich.
Optimization with PDE constraints, volume 23.
Springer Science & Business Media, 2008.

blauth2021optimal
Sebastian Blauth, Christian Leithäuser, and René Pinnau.
Optimal control of the sabatier process in microchannel reactors.
Journal of Engineering Mathematics, 128(1):1–28, 2021.

butt2012optimal
AIK Butt and R Pinnau.
Optimal control of a non-isothermal tube drawing process.
Journal of Engineering Mathematics, 76(1):1–17, 2012.

gotz2010optimal
Thomas Götz and SSN Perera.
Optimal control of melt-spinning processes.
Journal of Engineering Mathematics, 67(3):153–163, 2010.

ou2011unsteady
Kui Ou and Antony Jameson.
Unsteady adjoint method for the optimal control of advection and
  burger's equations using high-order spectral difference method.
In 49th AIAA Aerospace Sciences Meeting including the New
  Horizons Forum and Aerospace Exposition, page 24, 2011.

he2020lagrange
Ji-Huan He.
Lagrange crisis and generalized variational principle for 3d unsteady
  flow.
International Journal of Numerical Methods for Heat & Fluid
  Flow, 30(3):1189–1196, 2020.

wang2020variational
KangLe Wang.
Variational principle and its fractal approximate solution for
  fractal lane-emden equation.
International Journal of Numerical Methods for Heat & Fluid
  Flow, 31(7):2279–2287, 2020.

he2007variational
Ji-Huan He.
Variational iteration method—some recent results and new
  interpretations.
Journal of computational and applied mathematics, 207(1):3–17,
  2007.

berkani2012optimal
Samira Berkani, Farida Manseur, and Ahmed Maidi.
Optimal control based on the variational iteration method.
Computers & Mathematics with Applications, 64(4):604–610,
  2012.

akkouche2014optimal
Abderrahmane Akkouche, Ahmed Maidi, and Mohamed Aidene.
Optimal control of partial differential equations based on the
  variational iteration method.
Computers & Mathematics with Applications, 68(5):622–631,
  2014.

casas1995optimal
Eduardo Casas, Luis A Fernández, and Jiongmin Yong.
Optimal control of quasilinear parabolic equations.
Proceedings of the Royal Society of Edinburgh Section A:
  Mathematics, 125(3):545–565, 1995.

wang2009variational
Shu-Qiang Wang.
A variational approach to nonlinear two-point boundary value
  problems.
Computers & Mathematics with Applications,
  58(11-12):2452–2455, 2009.

he2023variational
Chun-Hui He and Chao Liu.
Variational principle for singular waves.
Chaos, Solitons & Fractals, 172:113566, 2023.

huang2022acceleration
Yakui Huang, Yu-Hong Dai, Xin-Wei Liu, and Hongchao Zhang.
On the acceleration of the barzilai–borwein method.
Computational Optimization and Applications, 81(3):717–740,
  2022.

zhang2012modelling
Xuan Zhang, Christian T Veje, Benny Lassen, and Morten Willatzen.
Modelling of the heating process in a thermal screw.
In Journal of Physics: Conference Series, volume 395, page
  012157. IOP Publishing, 2012.

huerre1985absolute
Patrick Huerre and Peter A Monkewitz.
Absolute and convective instabilities in free shear layers.
Journal of Fluid Mechanics, 159:151–168, 1985.




