Research article Institute of Brewing & Distilling

Received: 9 September 2013 Revised: 2 July 2014 Accepted: 3 July 2014 Published online in Wiley Online Library: 26 September 2014

(wileyonlinelibrary.com) DOI 10.1002/jib.150

Optimal control of beer fermentation processes

with Lipschitz-constraint on the control
Torsten Bosse* and Andreas Griewank
Nowadays, one can find in almost all industrial products a trail of mathematical optimization. In particular, the theory and
algorithms of optimal control have helped in various fields to reduce the production time and to improve the quality of
the considered products. As a special class of applications, two optimal control formulations for the fermentation process
of beer are presented. The reactions of the fermentation processes are modelled by a system of ordinary differential
equations that are steered by the heating and cooling of the involved substrates. The occurring physical limitations, such
as temperature limits and the requirement to avoid scenarios with unrealistic bang-bang controls, give rise to optimal control
problems with general control constraints. Therefore, this paper has reviewed a forward–backward sweeping method for
solving these kinds of optimal control problems and presents encouraging numerical results that were obtained by this
approach. Copyright © 2014 The Institute of Brewing & Distilling

Keywords: fermentation process of beer; optimal control with control constraints; forward–backward sweeping methods;
numerical results

Introduction variable functions y(t) = (y1(t), …, yn(t)), n ∈ , which satisfy a

given system of ordinary differential equations
Mathematical optimization is often a long and tedious task. It
only becomes interesting after constructing a good model for dy ðt Þ
which can one prove the existence of a solution and which has ¼ ψ ðy ðt Þ; uðt Þ; t Þ with initial values y ðt 0 Þ
dt
a solver at hand that can find one of those optimal solutions. ¼ y 0 ∈ℝn (ODE)
In particular it becomes interesting if the model represents
something that is beloved by mathematicians, such as coffee over the time interval I = [t0, tf] ⊂ ℝ and minimize the objective
or beer. In this paper we focus on the latter and review two functional
formulations for the fermentation process of beer (1–5) in terms Z tf
of optimal control problems that are steered by the heating and J ðy; u; t 0 ; tf Þ ¼ lðy ðt Þ; uðt Þ; t Þdt þ ϕ ðy ðtf Þ; uðt f Þtf Þ∈ℝ:
cooling of the involved substrates. The problems are solved by t0

the algorithm given in the literature (6,7) yielding an improved

solution that is applicable from a practical point of view. To simplify the subsequent discussion, we will avoid any func-
Moreover, two fermentation processes can be regarded as tion space analysis and assume the underlying functions to be
representative of almost all such applications. Hence, the continuously differentiable with regard to all variables over the
proposed algorithm can be used to implement a proportional- time interval I = [t0, tf]:
integral-derivative (PID) controller for small to large scale brew-
eries, for example, on a Raspberry PI computer (8). ψ ∈ C1;1 ðℝn ℝm I; ℝn Þ; ϕ∈C1;1 ðℝn ℝm I; ℝÞ;
We will consider a generalization of the unconstrained optimal
control problem: and l∈C1;1 ðℝn ℝm I; ℝÞ:
Z tf
min l ðy ðtÞ; uðt Þ; t Þdt þ ϕ ðy ðt f Þ; uðtf Þ; tf Þs:t: A more precise description can be found in the literature
y;u t0
(UCP) (9–11). Furthermore, we say that the optimal control problem
dy ðt Þ is unconstrained if the feasible set U(t) is equal to ℝm for all
¼ ψ ðy ðt Þ; uðt Þ; t Þ for all t∈½t 0 ; t f  and y ðt 0 Þ ¼ y 0
dt t ∈ I and constrained, otherwise.
By standard arguments from variational calculus one finds
which includes the additional control constraint
that any optimal solution (if existing – the existence can be
assured by additional properties on the underlying functions
uðtÞ∈Uðt Þ for all t∈½t 0 ; tf : (CC)

* Correspondence to: Torsten Bosse, Institut für Mathematik, Humboldt

Here the feasible set of the controls U(t) is assumed to be a Universität zu Berlin, Unter den Linden 6, 10099 Berlin, Germany. E-mail:
bosse@math.hu-berlin.de
nonempty convex set, possibly varying with the time t ∈ I. In
detail, we are looking for a vector of control functions u Institut für Mathematik, Humboldt Universität zu Berlin, Unter den Linden
(t) = (u1(t), …, um(t)) ∈ U(t) ⊆ ℝm, m ∈ , and a vector of state
444

6, 10099 Berlin, Germany

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

such as convexity or likewise) (y*, u*) of the unconstrained opti- Next, one applies a backward sweep to determine the
mal control problem UCP must be a stationary point, that is, it corresponding solution λ of the adjoint equation
satisfies the following stationarity conditions that are given in
terms of the corresponding Hamiltonian function
dλΤ ðt Þ ∂H0
¼ ðy ðtÞ; uðt ÞλðtÞ; tÞ with λðtf Þ
dt ∂y
H0 ðy ðt Þ; uðt Þ; λðt Þ; t Þ ¼ l ðy ðtÞ; uðtÞ; tÞ (HAM0) ∂ϕ
¼ ðy ðtf Þ; uðt f Þ; λðt f Þ; tf Þ:
þλΤ ðt Þψ ðy ðtÞ; uðtÞ; tÞ: ∂y

by integrating backward in time over I from tf to t0. In the last

Theorem (9) Let u* be a vector of optimal control functions and step, one then adapts the control for u(t) each time t ∈ I in such
y* be a vector of corresponding optimal state variables for the a way that the objective functional is reduced, based on the
unconstrained problem UCP. Then there exists a piecewise
  solution of the original system and on the computed sensitivity
differentiable vector-valued function λ ðtÞ ¼ λ1 ðt Þ; …; λm ðt Þ information from forward sweep and the backward sweep,
such that for all t ∈ I respectively. These three steps are repeated in an outer iteration
until the state, adjoint, and control variables satisfy the neces-
∂H0 sary optimality conditions FONC up to a certain precision or
0¼ ðy ðt Þ; uðt Þ; λðtÞ; tÞ where some other stopping criteria are satisfied, for example,
∂u
the  2 -norm of the control update Δu(t) over the time interval
dy  ðtÞ ∂H0 
¼ ðy ðtÞ; u ðtÞ; λ ðtÞ; tÞ ¼ ψ ðy  ðt Þ; u ðt Þ; t Þ; (FONC) I is smaller than some prescribed tolerance or a maximum num-
dt ∂λ
ber of iterations is exceeded.
dλ ðt Þ ∂H0 
¼ ðy ðt Þ; u ðt Þ; λ ðtÞ; tÞ; Using this approach, Gee and Ramirez (5) found an optimal
dt ∂y solution (bang-bang) for the second fermentation problem
(presented below) with box-constraints on the control, that is,
with initial values y  ðt 0 Þ ¼ y 0 and λ ðtf Þ ¼ ∂ϕ  
∂y ðy ðt f Þ; u ðt f Þ; t f Þ.
for each t ∈ I the control u(t) was limited from below and above
Assuming that an optimal solution exists, there is a huge by some constant. A similar approach was pursued by Griewank
variety of numerical methods to solve the problem UCP or at and Radwan (6) and Radwan (7) for the more general optimal
least find some stationary point. A common class of methods problem with control constraints
are the so-called sweeping methods (6,7,9,11) that are based on Z tf
Pontryagin’s maximum principle (12). The most simple sweeping miny;u lðy ðtÞ; uðt Þ; t Þ dt þ ϕ ðy ðt f Þ; uðtf Þ; tf Þ s:t:
method is the forward–backward algorithm (9), which consists of t0
8
three parts that are applied consecutively in each iteration, < dy ðtÞ ¼ ψ ðy ðtÞ; uðtÞ; tÞ:
>
(OCP)
namely a forward sweep, a backward sweep and a control update. for all t∈½t0 ; t f  dt
In the forward sweep one computes a solution of the underlying >
:
uðt Þ∈Uðt Þ
system of ODEs (ODE) by some standard integration scheme
over the time interval I and a current control u. y ðt 0 Þ ¼ y 0

Table 1. First-order method for optimal control problems with control constraints (6,7)

1 Initialize feasible control uþ ðt Þ ¼ uðt Þ∀t∈½t 0 ; t f ,

2 β0 ∈ ℝ+, ρ ∈ (0, 1), TOL, MAXITER, J = ∞, γ, ITER=0
3 DO
4 IF (ITER>0)
5 Solve LocP with β for all y(t), u(t), λ(t), t∈ [t0, tf] → Δu(t)
6 Set u+(t) = u(t) + Δu(t)
7 ENDIF
8 Integrate forward from t0 to tf:
dy þ  
9 dt ðt Þ ¼ ψ y þ ðt Þ; uþ ðt Þ; t y þ ðt 0 Þ ¼ y ðZ
t 0 Þ with
  tf  
10 Evaluate J þ ¼ ϕ y þ ðtf Þ; uþ ðtf Þ; tf þ l y þ ðt Þ; uþ ðt Þ; t dt
t0
11 IF(J+ > J) THEN (Reject Update)
12 Increase β ¼ βρ
13 ELSE (Accept Update)
14 Set u(t) = u+(t), y(t) = y+(t), J = J+ and β = ρβ
15 Integrate backward from tf to t0:
dλΤ ðtÞ
16 dt ¼  ∂H0 ðyðtÞ;u∂yðtÞ;λðtÞ;tÞ with initial values λðt f Þ ¼ ∂φðyðtf Þ;u∂ðtyf Þ;λðtf Þ;tf Þ
17 ENDIF
18 ITER=ITER+1
19 WHILE ((ITER < MAXITER)and(‖Δu‖ > TOL))
20 Output y(t), u(t), J, …
445

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

Therefore, the authors proposed a forward–backward The production model itself is governed by several ordinary
sweeping method to solve the problem OCP instead of solving differential equations that simulate the fermentation process.
a large nonlinear problem of the discretized version of OCP or The first model that we consider is the one from De Andres-Toro
applying some genetic optimization algorithm as was done in et al. (1) and (2) given by the system:
Ramirez and Maciejowski (3) and de Andres-Toro et al. (1,13),
respectively. Hence, the algorithm presented in Table 1 can
Lag Phase y active þ y lat ¼ constant¼0:48y init
be regarded as a generalization for the case of optimal control
problems with control constraints of the method used in Gee. dy lat dx active ðt Þ
¼ ¼ μlag y lat
and Ramirez (5). dt dt
The main ingredient of the algorithm is to compute at each dy active
iteration, that is, after one forward and backward sweep, an Fermentation ¼ μy y active  k m y active þ μlag y lat
dt
update Δu(t) that defines a new control u+(t) = u(t) + Δu(t) for dy bottom
the next iteration whenever there is an improvement in the ¼ k m y active  μD y bottom
dt
objective functional:
ds
¼ μs y active
  dt
J þ :¼ J y þ ; uþ ; t 0 ; t f < J :¼ J ðy; u; t 0 ; t f Þ: de
¼ μa f y active ;
dt
μy 0 s 0:5sinit μD0
The control updates are computed by solving for each time where μy ¼ μD ¼
discretization point t ∈ [t0, tf] a small, local nonlinear problem 0:5sinit þ e 0:5sinit þ e
depending on the current state y(t), the control u(t), the adjoint e μ 0s μ 0s
f ¼1 μ ¼ s μ ¼ a :
λ(t) and some non-negative parameter β ∈ ℝ: 0:5sinit s k s þ s a k s þ s
(Ferm 1a)
Δuðt Þ ¼ arg minΔuðtÞ H0 ðy ðt Þ; uðt Þ þ ΔuðtÞ; λðtÞ; tÞ þ βkΔuðtÞk2
subject to uðt Þ þ Δuðt Þ∈U:
These equations and the parameters, which we will present
(LocP) below, describe the whole process that can be controlled by
the temperature T = T(t) over the time t ∈ [t0, tf], that is, u(t) = T(t)
In fact the parameter will be adapted throughout the in OCP. As is common, the process is divided into two consec-
algorithm and has two effects for our method. First, it guaran- utive parts: a lag and a fermentation phase. Within the fermenta-
tees the existence of solutions for the local nonlinear problems tion phase one distinguishes between three different types of
if β is chosen sufficiently large. Second, it acts as a penalty yeast cells, lag, active and dead, whose concentrations are de-
parameter that penalizes large control steps and, therefore, noted by the variables ylat, yactive and ybottom (see also Table 2).
can be regarded as a step size regularization similar to overesti- Based on its amount, the active yeast cells then transform the
mation methods (compare with Conn et al. (14). Starting from given sugar s into ethanol e. Within this process there emerge
some initial β0 ≥ 0, the parameter β gets increased in every so called by-products, which contribute odours and flavours.
iteration by a certain factor, for example, doubled (ρ = 0.5), if The latter two properties are mainly due to the emerging
the control update direction Δu does not yield an improvement diacetyl diac and ethyl acetate acet, which both follow the
in the objective function for the corresponding state and adjoint system of ordinary equations
variables. In some cases it is beneficial to allow β to be reduced
to avoid the method from making too slow progress owing to a
too high penalization on the step length. For example, one could Table 2. Nomenclature of the first fermentation process
reduce β after all successful iterations, where the objective
function was significantly reduced. However, the method is not Parameter Description Unit
guaranteed to converge to an optimal solution if one allows a
yactive Suspended active biomass g/L
reduction of β.
ylat Suspended latent biomass g/L
The rest of the paper is structured as follows: first we review the
yinit Initial suspended biomass g/L
two fermentation models from the literature (1–5,15). Then we
ybottom Suspended dead biomass g/L
modify the given models by including some additional Lipschitz- sinit Initial sugar (=130) g/L
continuity constraints to avoid unrealistic control scenarios like s Concentration of sugar g/L
bang-bang controls and give the necessary adaptations for the e Ethanol concentration g/L
algorithm. All the results obtained by our method for both acet Ethyl acetate concentration ppm
fermentation processes are presented in the final section. diac Diacetyl concentration ppm
μy Specific rate of growth
Fermentation process I μD Specific settle down rate
μs Specific substrate consumption
As described in the literature (1–5), the goal for the production μa Specific rate of ethanol production
of beer is twofold. The main objective is to obtain a product with f Fermentation inhibition factor
an excellent flavour at the end. On the other hand, the producer kdc Appearance rate
wants to have a large profit, that is, the production time and kdm Reduction or disappearance rate
arising costs should be minimal.
446

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

By-products dðacetÞ ds The resulting optimal control problem over the time interval
¼ μeas ¼ μeas μs y active
dt dt [t0, tf] = [0, 160] is of Mayer form
(Ferm 1b)
dðdiacÞ
¼ k dc sy active  k dm ðdiacÞe minu ϕ ðy ðtf Þ; uðt f Þ; tf Þ subject to
dt 8
< dy ðt Þ ¼ ψ ðy ðt Þ; uðt Þ; t Þ
>
The model parameters are given in terms of the Arrhenius func-
for all t∈½t0 ; tf  dt (1)
tion that reflects the temperature dependence of the reaction rate >
:
constant, and therefore, rate of a chemical reaction itself: uðtÞ∈U
y ðt0 Þ ¼ y ðt 0 Þ
31934:09 38313
μy 0 ¼ e 
108:31
km ¼ e 130:16

T þ 273:15 T þ 273:15 where u(t) denotes the control variable in a feasible domain
26589 10033:28 U, that is, the temperature T (in Celsius) must be within
μeas ¼ e89:92  μD0 ¼ e33:82 
T þ 273:15 T þ 273:15 U = [∞, 15] for all times t ∈ [0, 160]. In fact, we can shrink
11654:64 1267:24 the set of feasible controls to some smaller interval to avoid unre-
μs0 ¼ e41:92 þ μa0 ¼ e3:27 
T þ 273:15 T þ 273:15 alistic scenarios, for example, U = [0, 15] owing to the freezing
9501:54 34203:95 point of water. The objective functional is given by the modified
μlag ¼ e30:72  k s ¼ e119:63 þ objective ϕ(y(tf), u(tf), tf) = J1  3 = J1 + J2 + J3 and ψ(y(t), u(t), t) repre-
T þ 273:15 T þ 273:15
k dc ¼ 0:000127672 k dm ¼ 0:00113864 sents the formulae (Ferm 1a), (Ferm 1b) and (Ferm 1c) of the fer-
(Ferm 1c) mentation process with initial values y0 = y(t0) at starting time t0.

As an objective function, De Andres-Toro et al. (1) minimized

the weighted sum J1  4 = J1 + J2 + J3 + J4 of four terms that Fermentation process II
characterize a full-bodied and great tasting beer: Another model for the fermentation of beer was given by Gee
and Ramirez (5). It was based upon previous works by Ramirez
Measure the final ethanol production J 1 ¼ 10: and Maciejowski (3), Gee and Ramirez (4) and Engasser et al.
Measure the final diacetyl production J 2 ¼ þ5:73e95ðdiacÞ11:51 (16). The model consists of three parts, being the growth model,
Limit level of ethyl acetate at the end J 3 ¼ þ1:16e46ðacetÞ66:77 a nutrient and a flavour model that are given by the formulae
Z tf (Ferm 2a), (Ferm 2b) and (Ferm 2c), respectively:
Temperature limit along the process J4 ¼ þ μLB dt
t0 8
> dX
(Ferm 1d) >
> Biomass ¼ μx X
>
> dt
>
>
>
> dG
>
< Glucose ¼ μ1 X
It was not exactly stated in the above mentioned paper what dt
Growth (Ferm 2a)
the integrand of the last formula was. From the text we >
> dM
>
> Maltose ¼ μ2 X
concluded that it was given by μLB = max(0, T  15) (in Celsius), >
> dt
>
>
which penalizes a too high temperature over 15 °C, which >
>
: Maltotriose dN ¼ μ X
increases the spoilage risk of the fermentation results. In an dt 3
unpublished communication (Gabrielle E. Carrillo with P. D.
Roberts and V. M. Becerra, Optimal control of fermentation 8
>
> dL dX L
processes, 1999, personal communication) that first brought our >
> Leucine ¼ Y LX D
>
> dt dt K L þ L
attention to this topic, the authors aimed to reduce the >
< dI dX I
production time to 160 h with a similar or better objective value Nutrient Isoleucine ¼ Y IX D (Ferm 2b)
>
> dt dt K I þ I
starting from a given temperature profile over 200 h with >
>
objective value Jinit ≈ 533.9. Therefore, they optimized the >
> dV dX V
>
: Valine ¼ Y VX D
temperature of the process using three different algorithms dt dt K V þ V
yielding an optimal profile with an objective value Jopt ≈ 568.3.
However, the resulting ‘optimal’ temperature profiles (1) were 8
> dIB dIA
>
> ¼ Y IB μV X; ¼ Y IA μL X
more than disappointing since they were not practical, that is, >
> dt dt
>
>
the optimal control had several jumps that cannot be achieved >
>
>
> dMB dP
in a real production situation. Moreover, the profile exhibited >
> ¼ Y MB μI X; ¼ Y PE ðμV þ μI ÞX
>
> dt dt
several kinks that surpassed the critical temperature of 15 °C, >
>
>
>
which might result in a spoiled end product. The same observa- >
> dEA
>
< ¼ Y EA ðμ1 þ μ2 þ μ3 ÞX
tions held true for the results in De Andres-Toro et al. (1), where dt
Flavour
the authors used a genetic algorithm to compute an optimal >
> dEC dIAc
>
> ¼ Y EC μX X; ¼ Y IAc μIA X
temperature profile of 130 h. >
>
>
> dt dt
For this reason, we now apply our optimal control solver to the >
>
>
> dVDK
same problem but with the following minor modifications. Instead >
> ¼ Y VDK μX X  k VDK ½VDK X
>
>
of considering the original objective function (Ferm 1d), we now >
> dt
>
>
minimize J1  3 = J1 + J2 + J3 and include the additional constraint >
>
: dAAl ¼ Y AAl ðμ1 þ μ2 þ μ3 ÞX  k AAl ½AAl X
T ≤ 15 °C to guarantee that the temperature always stays below dt
the critical value. (Ferm 2c)
447

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

 
E μi
dT 1 μi ¼ μi0 exp  ; i ¼ G; M; N
¼ ððμ1 ΔHFG þ μ2 ΔHFM þ μ3 ΔHFN ÞX þ pðT c  T ÞÞ RðT þ 273:15Þ
dt ρcp  
E Ki
(Ferm 2d) K i ¼ K i0 exp  ; i ¼ G; M; N; L; V (Ferm 2f)
RðT þ 273:15Þ
!
E ′K i
Similar to the first model, a lag effect described by the term K ′i ¼ K ′i0 exp  ; i ¼ G; M
RðT þ 273:15Þ
D ¼ 1  et=τD was included for the nutrients to simulate the
time delay that is possibly due to the time needed to transport
the amino acids across the cellular membranes of the yeast. Also The corresponding parameters of the system above are given
three classes of components that mainly influence the taste of by the following equations and constants listed in Table 3.
the fermentation result are modelled within the flavour model.
The first class contains the undesirable fusel alcohols, namely μG G μ M K ′G μ N K ′G K ′M
μ1 ¼ ; μ2 ¼ M ′
; μ3 ¼ N ′ ′
;
isobutyl alcohol (IB), isoamyl alcohol (IA), 2-methyl-l-butanol KG þ G K M þ MK G þ G K N þ N K G þ GK M þ M
(MB), and propanol (P). The good flavour components form the KX
second class are the esters: ethyl acetate (EA), ethyl caproate μX ¼ ðμ1 Y XG þ μ2 Y XM þ μ3 Y XN Þ ;
K X þ ðX  X 0 Þ2
(EC) and isoamyl acetate (IAc). In the last class we find the vicinal
1 dV 1 dL 1 dI
diketones (VDK) and acetaldehyde (AAl) that cause a bad flavour and μV ¼  ; μL ¼  ; μI ¼  :
X dt X dt X dt
and, therefore, are also considered as undesirable quantities.
Moreover, an additional equation was included to mimic the tem-
perature of the process emerging from the chemical reactions. The Most of the values in Table 3 coincide with the original ones
control is now given by a user supplied cooling function Tc, that is, given in Gee and Ramirez (5). However, we were not able to find
Tc(t) = u(t) regulating the current temperature T(t) via (Ferm 2d), all needed quantities within the cited paper and, therefore, were
where p ¼ pp1 p2 with p1 denoting the constant heat transfer coeffi- forced to take appropriate values from additional literature (3,4),
3
cient (kJ/hm2C∘), p2 being the heat transfer area (m2), and p3 being which might result in inconsistencies to the referred literature.
the volume of the fermenting mixture (m3). Further quantities such Furthermore, we needed to take a guess for the initial values
as the amount of produced ethanol are given by the equation of the ODE from the figures in Gee and Ramirez (5):

∂E y ðt0 Þ ¼ ½X 0 ; G0 ; M0 ; N0 ; L0 ; I0 ; V 0 ; IB0 ; IA0 ; MB0 …

¼ X ðY EG μ1 þ Y EM μ2 þ Y EN μ3 Þ: (Ferm 2e) ¼ ½125; 70; 220; 40; 1:2; 0:6; 1; 0; 0; 0…
∂t
and the specific growth rates depending on the temperature
being of Arrhenius-type analogous to the first fermentation As the question mark in Table 3 indicates, we also did not find
process, appropriate values for the estimation of parameter μIA required

Table 3. Parameter estimates of the second fermentation process according to Ramirez and Maciejowski (3)

Parameter Value Parameter Value Parameter Value

KX 365000 X0 125
G0 70 M0 220 N0 40
ln μG0 35.77 ln μM0 16.4 ln μN0 10.59
EμG 22.6 × 103 EμM 11.3 × 103 EμN 7.16 × 103
ln KG0 121.3 ln KM0 19.5 ln KN0 26.78
EKG 68.6 × 103 EKM 14.4 × 103 EKN 19.9 × 103
ln KL0 10.14 ln KV0 328.0 KI 0.365
EKL 5.95 × 103 EKV 211.9 × 103
lnK ′G0 23.33 lnK ′M0 55.61
E ′KG 10.2 × 103 E ′KM 26.3 × 103
YXG 0.134 YXM 0.268 YXN 0.402
YEG 1.92 YEM 3.84 YEN 5.76
YLX 0.0832 YIX 0.0363 YVX 0.0273
YIB 0.203 YIA 0.557
YMB 0.472 YPE 0.235
YEA 0.000992 YEC 0.000118 YIAc 0.0269
YVDK 0.000105 YAAl 0.01
lnk 0VDK 86.8 lnk 0AAl 10.4
EVDK 54.3 × 103 EAAl 11.1 × 103
ΔHFG 91.2 ΔHFM 226.3 ΔHFN 361.3
ρ 1040 cp 4.016 R 1.987
τD 23.54 μIA ???
p1 1.6 p2 0.188 p3 0.1
448

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

for the equation describing the behaviour of the isoamyl acetate characteristics of the vessel, the heat exchanger unit and the
(IAc). Thus, we decided to remove the corresponding formula refrigeration load during the heat transfer pose constraints on the
from the model since it is not important for the optimization temperature profiles that can be realized during the fermentation
procedure, that is, neither of the other equations was affected process. In particular, one expects that the heating/cooling rate of
by this quantity nor the considered objective function: the brine temperature is limited by some parameter γ and should
be taken care of in the optimization of the process.
J ¼ E ðt f Þ þ C 1 ½IBðtf Þ þ IAðtf Þ þ MBðt f Þ þ Pðt f Þ (Ferm 2h) Therefore, we consider a modification of the control con-
straint in the following. For fixed γ ∈ ℝ+, we define the modified
þC 2 AAl ðtf Þ
e
feasible set of controls U⊆U to be
where C1 = 810 and C2 = 610.  
Ue ¼ u∈C ð½t0 ; tf Þ with Lipschitz constant γ and u∈U (Lip)
0;1
As in the first model we had a problem of Mayer form with
objective J ¼ ϕ tf (Ferm 2h) and ψ representing (Ferm 2a)– that only allows temperature profiles in the original feasible set
(Ferm 2d) without the equation for isoamyl acetate (IAc) over U, for example, in the interval [Tmin, Tmax], with a maximal ascent
the time interval [t0, tf] = [0, 150]. or minimal descent rate γ. To realize this requirement, we
included the extra constraint

Lipschitz constraints on the control juðtiþ1 Þ  uðti Þj≤ γjt iþ1  t i j

Within a fermentation process the temperature of the brine is to the local problem, where ti < ti + 1 denotes two consecutive
usually controlled by the brine jackets attached to the conical time-steps ti, ti + 1 ∈ [t0, tf] from the discretization scheme. Thus,
vessel that contains the fermentation mesh. The physical if u(ti) was the previous control with corresponding update Δu(ti)

Figure 1. Reference temperature profile for example 5.1 and corresponding state variables over a time interval of 160 h. This figure is available in colour online at
449

wileyonlinelibrary.com/journal/brewing.

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

at time ti then the local problem LocP for the calculation of the Note that, for the second fermentation process, one must
control step Δu(ti + 1) at time ti + 1 reads ensure that there always exists a solution, for example, a too
min H0 ðy ðt iþ1 Þ; uðtiþ1 Þ þ Δuðt iþ1 Þ; λðtiþ1 Þ; tiþ1 Þ þ βkΔuðt iþ1 Þk2 restrictive choice of the cooling rate might cause a heating effect
Δuðtiþ1 Þ
of the chemical reaction that cannot be compensated for by the
s:t: T ðt iþ1 Þ þ Δuðtiþ1 Þ in Uðtiþ1 Þ; cooling mechanism, and the temperature of the whole process
and Δuðtiþ1 Þ≤uðti Þ þ Δuðt i Þ þ γðt iþ1  t i Þ  uðtiþ1 Þ could increase too rapidly.
Δuðt iþ1 Þ≥uðt i Þ þ Δuðti Þ  γðtiþ1  ti Þ  uðtiþ1 Þ
Numerical results
In other words, the Lipschitz constraint on the control can be For all examples, we used Euler’s explicit one-step method
regarded as a varying box-constraint on Δu(ti + 1) depending on to simulate the fermentation process on an equidistant time
the computed control u+ (ti) = u(ti) + Δu(ti) at the previous time discretization with time step length Δt = 2 h. Of course, we
step ti. Note that the above formulations fit perfectly the way could also have used higher integration methods, adaptive
the considered algorithm is designed, where the steps Δu(ti) discretization schemes and a finer grid but in our testing we
are computed consecutively. As can be seen in the following found that this choice was already sufficient to give a reason-
section containing the numerical results, the proposed able result. The code was implemented in MATLAB using
algorithm finds an optimal solution yielding a control profile that FMINCON’s active-set strategy for the local sub-problems and
behaves in a friendlier manner, that is, there are not that many ADIMAT (17) for the computation of the necessary derivatives.
jumps as in Ramirez and Maciejowski (3), where some smooth- (A current version of the code can be obtained upon request
ing techniques are required. However, the result is still a bang- from the corresponding author. The software is free of any
bang control that is hard to realize from a practical point of view. guarantee and is only allowed for academic use.) The

Figure 2. Solution for example 5.1 in the case 0 ≤ u(t) ≤ 15 and u′(t) ∈ ℝ. This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.
450

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

tolerance of FMINCON was set to 108, whereas we used a Note that this updating procedure differs from the proposed
maximal number of iterations of ≤ 1000, a maximal β < 1016 strategy, where the authors suggested choosing β = 0 in the first
and a minimal step length kΔ uk 2 > 1016 as stopping step and then starting to increase β if no reduction in the objective
criteria for the optimal control method. The factor β was functional is found. However, in our numerical tests the update of β
initialized with β0 = 1 and adapted according to the following (given above) proved to be efficient and still reliable. The initial con-
strategy: trols for both problems were chosen to be the reference profiles.

(  
max 1016 ; β=2 if reduction in the previous two iterations
βnew ¼ :
β ¼ 10β if no reduction in the previous iteration

Table 4. Results for the reference temperature profile (Fig. 1), the box-constrained optimization (example 1.1) with γ = 0.1
(example 1.2) and γ = 0.01 (example 1.3)

Example Bound γ J 13

opt J 1opt J 2opt J 3opt It
Reference / — 4.72 × 102 4.72 × 102 2.31 × 104 1.11 × 1028 1
1.1 0/15 — 5.98 × 102 5.98 × 102 7.31 × 105 4.19 × 101 153
1.2 0/15 0.1 5.91 × 102 5.91 × 102 1.78 × 104 1.96 × 101 14
1.3 0/15 0.01 4.99 × 102 4.99 × 102 2.02 × 104 3.97 × 1028 108

Figure 3. Solution for example 5.1 in the box-constrained case 0 ≤ u(t) ≤ 15 with |u′(t)| < γ = 0.1. This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.
451

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

Fermentation process I time interval I = [0, 160]. Here, the term bounds denotes the
For the first model, we used the temperature profile depicted in lower/upper bound on the control u(t) and ‘It’ the number of
Fig. 1 as reference control with objective value J1  3 ≈  472 over required iterations. From the values of the fourth column (J 13opt ),
the time interval I = [0, 160]. The optimal temperature profile for one can see that the optimization improves the reference profile
the box-constrained problem and the resulting quantities that in all three cases and that the less restrictive case of only having
were obtained by the proposed method are visualized in Fig. 2. box-constraints yields better results than the more restrictive
The optimal function values can be found in Table 4. Note that problems, including the additional Lipschitz-constraints γ = 0.1
the obtained optimal control is feasible and yields a significantly and γ = 0.01, respectively.
improved objective value compared to the reference profile. All optimal profiles exhibit the expected behaviour, that is,
Moreover, the results are consistent in that the more restrictive at first a maximal heating is suggested followed by a cooling
problems do not allow for better solutions. phase, which is consistent with the results derived by the
However, as can be seen in Fig. 2 it might happen that the other authors cited. In particular, the box-constrained case
profiles are not practical owing to too large temperature depicted in Fig. 2 converges towards a bang-bang control,
variations that cannot be realized in practice. Therefore, we where the temperature jumps from approximately 14.5 to
considered the Lipschitz-constrained case as discussed above 10 °C at t ≈ 100. For the Lipschitz-constrained case, it can be
for the two arbitrary choices γ = 0.1 and γ = 0.01 that yield the seen (see Fig. 5) that the maximal rate γ of increase and
optimal solutions given in Figs 3 and 4, respectively. The decrease is attained by the optimal temperature profile at
corresponding ascend and descend rates of the control for all the beginning and the end of the process, respectively.
cases are depicted in Fig. 5. Furthermore, it can be observed that a too restrictive choice
A short summary of the computed numerical values is of γ results in a ‘heating–cooling’ strategy (see Fig. 4) with
presented in Table 4 for N = 65 discretization points over the an optimal switching point at t ≈ 125.

Figure 4. Solution for example 5.1 in the box-constrained case 0 ≤ u(t) ≤ 15 with |u′(t)| < γ = 0.01. This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.
452

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

Figure 5. Lipschitz constant for example 5.1 of the reference temperature (Ref) (top-left), the box-constrained control example 1 (top-right) and the two Lipschitz cases
γ = 0.1 (bottom-left) and γ = 0.01 (bottom-right). This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.

In contrast, larger values for γ yield a stationary phase at option since negative quantities do not have a physical meaning
constant temperature (see Fig. 3, 50 ≲ t ≲ 120). Moreover, the in this application and the results approve this modification.
constant temperature at the end of the process (see Fig. 3, The feasible set U for the control steps Δu(t) in the local nonlinear
120 ≲ t ≲ 160) suggests that for fixed parameter γ less time is LocP was defined by the box-constraints  u(ti) ≤ Δu(ti) ≤ 80  u(ti)
required to produce similar fermentation result. to limit the cooling rate TC(t) = u(t) to the interval [0, 80]. Moreover,
we added a penalty term to the objective
Z tf
Fermentation process II ϵ maxð0; T ðt Þ  15Þ2 dt ⇒ lðy ðt Þ; uðt Þ; t Þ ¼ maxð0; T ðt Þ  15Þ2
t0
In the second model, the cooling of the process is considered as
(Pen)
a control variable instead of the temperature profile itself. Here,
a reference profile was chosen similar to the first problem (cf. with ϵ = 10 to avoid temperatures T over 15 °C, respectively.
6

Fig. 6) as a constant cooling with Tc(t) = 10 for all t ∈ [0, 150] Analogous to the numerical results of the first fermentation
and J ≈ 324 as objective value (without the penalty term). process in the previous subsection, we present all optimal temper-
To achieve this result, it was necessary to include the minor ature profiles for the box-constrained case and the two Lipschitz-
correction within the ODE solver by only allowing non-negative constrained cases with γ = 0.1 and γ = 0.01 in the Figs 7–9. The
concentrations of the considered components, that is, we set the corresponding numerical values are listed in Table 5 and are
quantities of the involved substrates to the maximum of zero again consistent in that the less restrictive problems allow for
and their original value after each integration step. Note that this better optimal results.
is inconsistent with the requirement that the functions are Except of the reference profile, all solutions exhibit the same
continuously differentiable. Nevertheless, we preferred to use this characteristic behaviour, namely, a monotonic increase in the
453

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

Figure 6. Reference temperature profile (Ref) for example 5.2 and corresponding state variables over a time interval of 150 h. This figure is available in colour online at
wileyonlinelibrary.com/journal/brewing.
454

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

Figure 7. Solution for example 5.2 with 0 ≤ u(t) ≤ 80, u′(t) ∈ ℝ, and penalty term Pen. This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.
455

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

Figure 8. Solution for example 5.2 with 0 ≤ u(t) ≤ 80, |u′(t)| ≤ 0.1, and penalty term Pen. This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.
456

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458
Optimal control of beer fermentation processes
Institute of Brewing & Distilling

Figure 9. Solution for example 5.2 with 0 ≤ u(t) ≤ 80, |u′(t)| ≤ 0.01, and penalty term Pen. This figure is available in colour online at wileyonlinelibrary.com/journal/brewing.
457

J. Inst. Brew. 2014; 120: 444–458 Copyright © 2014 The Institute of Brewing & Distilling wileyonlinelibrary.com/journal/jib
 T. Bosse and A. Griewank
Institute of Brewing & Distilling

Table 5. Results for the reference temperature profile in Fig. 6, the box-constrained case (example 2.1) with γ = 0.1 (example 2.2)
and γ = 0.01 (example2.3)
R
Example Bound γ Jopt E(tf) C1[…] C2[…] l(y, u, t)dt It
Reference / — 2.07 × 10 7
1.18 × 10 3
1.24 × 10 3
2.62 × 10 2
2.07 × 10 1
1
2.1 0/80 — 7.78 × 102 9.00 × 102 1.22 × 103 4.48 × 102 9.15 × 1011 77
2.2 0/80 0.1 8.29 × 102 8.62 × 102 1.23 × 103 4.64 × 102 1.42 × 1012 43
2.3 0/80 0.01 9.72 × 102 6.68 × 102 1.19 × 103 4.47 × 102 2.04 × 106 16

temperature to the maximal allowed value of 15 °C at the final an industrial fermentation case, Proceedings of the 36th IEEE Confer-
time tf = 150. The temperature for the reference profile itself ence on Decision and Control, Vol. 1, pp. 828–829.
2. De Andres-Toro, B., Giron-Sierra, J., Lopez-Orozco, J., Fernandez-
violates this upper limit since the cooling cannot compensate Conde, C., Peinado, J., and Garcia-Ochoa, F. (1998) A kinetic model
for the heat emerging from the chemical reaction. The resulting for beer production under industrial operational conditions, Math.
optimized cooling profiles show again the typical ‘heating– Comput. Simul. 48, 65–74.
cooling’ behaviour that was also observed for the first 3. Ramirez, W. F., and Maciejowski, J. (2007) Optimal beer fermentation,
J. Inst. Brew. 113(3), 325–333.
fermentation process.
4. Gee, D. A., and Ramirez, W. F. (1994) A flavour model for beer
fermentation, J. Inst. Brew. 100(5), 321–329.
5. Gee, D. A., and Ramirez, W. F. (1988) Optimal temperature control for
Conclusions batch beer fermentation, Biotechnol. Bioeng. 31(3), 224–234.
6. Griewank, A., and Radwan, A. (2011) First order method for optimal
In this paper, we reviewed the forward–backward sweeping control using parametric optimization, Proceedings of the World
method given in Griewank. and Radwan (6) and Radwan (7) for Congress on Engineering, Vol. I.
solving optimal control problems with control constraints. The 7. Radwan, A. (2012) Utilization of parametric programming and
algorithm was applied to the two beer fermentation processes evolutionary computing in optimal control. PhD thesis, Humboldt-
presented in the literature (1–5,15). For both examples, we found Universität zu Berlin.
8. BrewPi (2013) Homepage of BrewPi. Available from: http://www.
a significant improvement in the objective functional for the brewpi.com (accessed June 2014).
computed optimal temperature profile compared to given 9. Lenhart, S., and Workman, J. (2007) Optimal Control Applied to Biolog-
reference profiles, even those from the referenced literature. ical Models, Chapman and Hall/CRC Mathematical and Computa-
For future research, we propose to implement the triple- tional Biology, London.
10. Vinter, R. (2010) Optimal Control, Modern Birkhäuser Classics,
sweep method presented in Radwan (7) and compare the Springer, Berlin.
efficiency of both algorithms for the considered problems. 11. Gerdts, M. (2012) Optimal Control of ODEs and DAEs, De Gruyter,
Furthermore, an application of the approaches to a real fermen- Berlin.
tation process is envisioned in order to verify the methods and 12. Lebedev, L., and Cloud, M. (2003) The Calculus of Variations and
models themselves. Also, an automatic online monitoring tool Functional Analysis: With Optimal Control and Applications in
Mechanics, Series on Stability, Vibration and Control of Systems,
[compare (8)] can be realized since one can hot-start the Series A, Vol. 12, World Scientific, Singapore.
methods at any time ti ∈ [t0, tf] of the fermentation process, 13. de Andres-Toro, B., Giron-Sierra, J. M., Lopez-Orozco, J. A., and
where only the initial values need to be set as the current data Fernandez-Conde, C. (1997) Optimization of a batch fermentation
and the considered time interval needs to be adapted. There- process by genetic algorithms, Departamento de Informatica y
Automatica, Universidad Complutense de Madrid.
fore, a more efficient implementation in C/C++ or Fortran with 14. Conn, A., Gould, N., and Toint, P. (2000) Trust Region Methods,
ADOL-C and TAPENADE and check-pointing is recommended. MPS-SIAM Series on Optimization, Society for Industrial and Applied
Mathematics, Philadelphia, PA, USA.
Acknowledgements 15. Carillo-Ureta, G. E., and City University (2003) Optimal Control of
Fermentation Processes, City University, London.
We would like to thank the reviewers for their valuable comments 16. Engasser, J. M., Marc, I., Moll, M., and Duteurtre, B. J. M. (1981) Kinetic
and many helpful suggestions that allowed us to improve and modeling of beer fermentation, Proc. Eur. Brew. Conv. Congr., Copen-
clarify the manuscript. hagen, pp. 579–586, IRL Press, Oxford.
17. Bischof, C. H., Bücker, H. M., Lang, B., Rasch, A., and Vehreschild, A.
(2002) Combining source transformation and operator overloading
References techniques to compute derivatives for MATLAB Programs,
Proceedings of the Second IEEE International Workshop on Source
1. De Andres-Toro, B., Giron-Sierra, J., Lopez-Orozco, J., and Fernandez- Code Analysis and Manipulation (SCAM 2002), IEEE Computer Society,
Conde, C. (1997) Using genetic algorithms for dynamic optimization: Los Alamitos, CA, pp. 65–72.
458

wileyonlinelibrary.com/journal/jib Copyright © 2014 The Institute of Brewing & Distilling J. Inst. Brew. 2014; 120: 444–458