Using CMA-ES for Tuning Coupled PID Controllers
within Models of Combustion Engines

arXiv:1609.06741v4 [] 5 Jun 2017

Kateřina Henclová

1

Regional Innovation Centre for Electrical Engineering,
University of West Bohemia,
Univerzitnı́ 26, 306 14 Pilsen,
Czech Republic

Abstract
Proportional integral derivative (PID) controllers are important and widely
used tools in system control. Tuning of the controller gains is a laborious
task, especially for complex systems such as combustion engines. To minimize
the time of an engineer for tuning of the gains in a simulation software, we
propose to formulate a part of the problem as a black-box optimization task.
In this paper, we summarize the properties and practical limitations of
tuning of the gains in this particular application. We investigate the latest
methods of black-box optimization and conclude that the Covariance Matrix
Adaptation Evolution Strategy (CMA-ES) with bi-population restart strategy, elitist parent selection and active covariance matrix adaptation is best
suited for this task. Details of the algorithm’s experiment-based calibration
are explained as well as derivation of a suitable objective function. The
method’s performance is compared with that of PSO and SHADE. Finally,
its usability is verified on six models of real engines.
Keywords: PID controller, tuning coupled controllers, CMA-ES, black-box
optimization
1. Introduction
In a running combustion engine, one or more PID controllers ensure that
certain quantities (e.g. intake pressure or exhaust gas temperature) remain
1

Email address: henclova@karlin.mff.cuni.cz

Preprint submitted to Applied Soft Computing

June 7, 2017

constant or within given range. Since these quantities may be naturally
related and affect each other, controllers are often coupled and their gains
cannot be tuned independently of each other.
When engines are modeled, as in simulation software WAVE2 , controllers
may and need to be tuned using simulations with little or no knowledge about
the transfer function that describes the system. The complicated finiteelement model cannot be described by a simple formula and, moreover, it
comes with its own modeling and discretization errors.
Presently, manual work makes up a major part of the controller tuning
process. This lengthy procedure is based on trial and error and requires a
knowledgeable and experienced control engineer. For systems with a single
controller (or multiple but decoupled controllers), simple rules of thumb can
be employed (e.g. Ziegler-Nichols [43]). Similar, already-solved problems can
also provide a guideline. However, when having a complicated or unique system of coupled controllers, the complexity of the task makes it very difficult
to solve even for an experienced control engineer. Moreover, in our application of PID controllers within combustion engine models, other professionals
need to tune the controllers as well, creating the need for a simple-to-use,
robust tool.
Our goal is to formulate at least a part of the problem and deliver a
method that would eliminate or significantly lower the need for manual tuning. It should find a solution of a given problem within acceptable time
and with as little user interaction as possible. When combined with simple
tuning rules or educated guesses, our method is to use the provided solution
approximation as a starting point and quickly find a more refined solution.
Further testing of solutions found by our method (such as verifying their
robustness) is, however, still left for manual post-processing. In principle,
it is simple to replace the proposed objective function by its expected value
over parametric uncertainties. However, the uncertainties on the parametric
space are not available for the considered application and thus it is not a
part of this work. Moreover, such approach would significantly increase the
computational load.
The proposed method is used to design the controller gains for a predefined trajectory of the controller. For more complex system that requires the
2

WAVE is 1D engine and gas dynamics simulation software package developed by Ricardo Software [36] and used by the author of this paper.

2

use of gain scheduling [1] or fuzzy PID [33], the proposed method will be
used to design the gain in the operational point designed by the supervising engineer. The engineer will also be responsible for selection of the gain
scheduling or fuzzy rules.
The PID tuning problem with either one controller or multiple but decoupled or symmetric controllers can be and has been reformulated as a
black-box optimization problem and solved with an appropriate method.
Evolutionary algorithms have been successfully applied to many engineering problems [6] and they have also been used to tune PID controllers. The
genetic algorithm [10, 24] was one of the first to be used [30], and many others followed such as the differential evolution (DE) [35, 37, 38] in [5, 25] or
the particle swarm optimization (PSO) [28] and its combinations (e.g. with
bacterial foraging algorithm) in [8, 9, 25, 29].
The tuning problem with multiple coupled controllers can be formulated
as an optimization problem as well but it adds another level of complexity to
the task. (To the author’s knowledge, there have not been any papers dealing
with such problems.) Even though there are naturally multiple objectives
(one for each controlled quantity), we can use a simple trick to carefully
transform them into a single objective and thus enable the problem to be
solvable by usual means. This way, however, the problem becomes harder to
solve, while the time budget remains small. With simulations taking up to
several minutes each, we aim for an overnight or a one-day computation on a
regular PC, i.e. a few thousand simulation runs at most. This imposes high
expectations upon efficiency of the method used.
Considering properties of the problem, we choose to use a variant of the
Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [17, 21, 22,
27], an evolutionary algorithm founded deep in probability theory. It has
proven to be very effective and robust method in the extensive testing of
Black-Box Optimization Benchmarking (e.g. [20, 4]), surpassing the above
mentioned algorithms and many others (on the relevant sort of problems).
Despite its fame in the optimization community and large number of practical
applications, it has so far been little used for tuning PID controllers [25, 26,
42] or similar problems [19, 18].
In this paper, we experiment with the algorithm’s settings and come up
with a method that meets practical requirements, such as tolerable runtime,
and is thus fit for common use by engineers working with complicated simulations (e.g. users of WAVE). Moreover, we demonstrate its applicability on
six models of real-world engines with one, two and three controllers.
3

For comparison with CMA-ES, we use PSO (as implemented in the software package DEAP [7]) and a restarted version of SHADE (Success-History
based Adaptative DE [39]) with an implementation based on that by Tanabe
[40]. We show that they do not perform by far so well as CMA-ES.
2. Formulation of the problem
In this section, we formulate the controller tuning problem as a problem
of numerical optimization. This part is essential as it determines whether the
problem will be solvable (with the given practical limits) and what method
should be used. Compared to other research on controller tuning [5, 8, 9, 25,
29, 30], dealing with coupled controllers requires an extra level of complexity.
2.1. PID controllers in the context of engine simulations
PID controllers are well known and powerful tools in system control [11,
32]. To describe their inner workings, we will use the traditional, though
perhaps old-fashioned, notation. A controller’s input is the error function
e(t) = |actual(t) − target(t)| ,

(2.1)

i.e. the time-dependent absolute difference between the desired target value
and the actual value of a quantity (as measured by a sensor or computed
by a model). The output control signal that defines the system’s subsequent
reaction is given as
Z t
d
(2.2)
e(τ )dτ + D e(t),
C(t) = P e(t) + I
dt
0
where P , I and D are the proportional, integral and derivative gains, respectively.
The controllers’ implementation is provided within the simulation software WAVE. All the computations as well as the whole engine model are
given to us as a black box and we do not analyze them. The software as
well as the models are carefully calibrated to comply with the real world.
The simulation is fully deterministic and we do not introduce any random
perturbations of environment variables during the tuning process.
We input the controllers’ gains into the simulation software and obtain
computed trajectory of the controlled quantities over time. Our goal is to find
such constant gains P , I and D for each controller, that the corresponding
4

controlled quantities converge to the target value and do so as quickly as
possible. This provides us with the key how to compare solutions’ quality.
The entire tuning process is offline but, due to nonlinearities in the system, the curve is not as smooth as it is common for simple systems.
In this work, we are always provided with initial values and with (nonzero)
constant target values of the controlled quantities. The values are given by
the context of combustion engines.
2.2. The optimization problem
With k controllers within a system, each determined by three constant
gains P , I and D, there are 3k gains to be tuned: x = (P1 , I1 , D1 , . . . , Pk , Ik ,
Dk ). Let us further suppose that we are given an objective function M :
R3k → R. It returns a single number for every vector of controllers’ gains
and describes how well the whole system is controlled by controllers with such
gains. Assuming that higher quality inputs have lower function values, we
search for the objective function’s minimizer. The numerical optimization
problem derived from the controller tuning problem can be formulated as
looking for x∗ such that M (x∗ ) = essinf x∈R3k M (x). I.e. x∗ is a minimizer
of M that is reachable by numerical means. Naturally, we do not strive to
find the perfect optimizer x∗ but rather any good-enough solution x∗∗ for
which it holds M (x∗∗ ) − M (x∗ ) < , where  > 0 is a given tolerance. Such
“minimizer” x∗∗ is not unique but, for practical purposes, it need not be.
In our case, function M must incorporate a finite-element model of an
engine and process its results. The model’s input is the vector of controllers’
gains (a candidate solution) and it outputs development of the controlled
quantities’ values over time (the controllers are a part of the model). These
outputs then need to be processed in order to obtain a single number. Here,
it is essential to carefully encode information about the candidate solutions’
quality, so that their final values rank them well.
2.3. Definition of the objective function
When the controllers’ gains are set and the simulation is run, it outputs
the error functions ei (t) = ei (x, t), i = 1, . . . , n, described by (2.1). We
shall now focus on how to process ei so that the final function value contains
all information about quality of the input. In this section, the vector of
controllers’ gains x is arbitrary but fixed. So, for simplicity of notation, it is
left out from the functions’ arguments: we write ei (t) instead of ei (x, t) etc.

5

2.3.1. Controlling multiple objectives
Let us assume for a moment that once controllers’ gains are set, quality of
the setting is characterized by functions Fi (t) : R3k → R, i = 1, . . . , n, that
describe the n controlled quantities’ development (simulation time t being
given). That provides us with n objectives to be minimized. In our case,
the number of controllers is the same as the number of controlled quantities
(i.e. k = n) but, in general, this is not required. We expect k to be small,
typically k ≤ 5.
We shall combine these objectives into a single one. Using a weighted sum
(see e.g. [31]), we can define the
Pnobjective function describing the response of
the whole model as F (t) =
i=1 wi Fi (t), where wi is the weight constant
corresponding to the i-th controlled quantity. The purpose of constants wi is
to express our priorities and, more importantly, make all the objective functions Fi comparable, as they may significantly differ in range depending on
the corresponding units. Therefore, we set wi = pi /scalei , where pi describes
priority of the i-th objective and scalei contains knowledge of its typical value
or range. In our case, we set scalei = targeti > 0 to be the (constant and
nonzero) target value of the i-th controlled quantity. The targeti element is
the remainder of the integral under its (constant) curve, considering that all
the objectives are measured for the same length of time. In our experiments,
we shall always set priorities pi = 1, i = 1, . . . , n.
Thus we have obtained the way of combining multiple objectives into one:
F (t) =

n
X
i=1

1
Fi (t).
|targeti |

(2.3)

2.3.2. Step-response based objective function
It remains to define Fi , the function corresponding to the i-th objective
(controlled quantity), i = 1, . . . , n. The general framework is based on how
the step response looks like, see the idealized step response in figure 1.
The step response describes how the controlled quantity behaves after an
engine is started. The abrupt step change of the target value is the most
difficult part for the controllers to deal with. Once the engine is running, the
changes in the environment and the target values are not so sudden. Thus,
when we tune the controllers for the “worst case scenario”, there is a good
reason to believe (supported by empirical evidence) that they will be able to
control other cases as well. When the environment changes too much, e.g.

6

Figure 1: For a given setting of controllers, the error at time t is given as e(t) =
|actual(t) − target(t)|.

hitting the extremes of RPM (rotations per minute), the controller gains may
be changed using gain scheduling [1].
Functions Fi are naturally given by formula
Z t
Fi (t) =
T (τ ) E(ei (τ )) dτ, t0 ≥ 0,
(2.4)
t0

where t0 is the initial time (typically t0 = 0), T (τ ) is a function of time,
E(ei (τ )) is a function of error, ei (τ ) is the error as defined in (2.1). In our
case, t is given by the end of the simulation.
Function E(ei (τ )) describes our concern over the actual error size. On
the other hand, function T (τ ) characterizes how much we care about when
the particular error occurs. Its value at a given time can also be perceived as
weight assigned to the corresponding error function value. There are multiple
commonly used criteria of step response quality (IAE, ITAE, ITSE or ISE),
which all fit the framework described above, giving us common choices of
functions E(ei (τ )) and T (τ ). The well known ITAE criterion [12]
Z t
IT AE(t) =
τ |e(τ )| dτ,
(2.5)
0

fits our purpose the best, as we are highly concerned about when the error
happens: it is not so important in the beginning, while at the end, error
means nonconvergence.
This information must not get lost even after combining multiple objectives in (2.3). Therefore, to sort the candidate solutions better, we slightly
modify (2.5). First, we shift the time, so that the value of one second loses
its importance as a factor. Then, we can choose to penalize nonconvergent
7

solutions. We can choose to stress the importance of time by increasing the
degree of its polynomial from linear to quadratic. Or, instead, we can leave
out the beginning of the time interval as the large error that can be expected
there is – from the optimization viewpoint – only noise. Various possibilities
have been experimented with and the best option found was to set
Z t
(τ + 1)|ei (τ )| dτ t0 ≥ 0.
(2.6)
Fi (t) =
t0

Here, the beginning of integration may be shifted from 0 to arbitrary t0 .
Selection of the shift t0 along with its justification and influence upon performance is discussed later. The simulation length t is set appropriately and
is same for all controlled quantities.
2.3.3. The objective function
In previous steps, we have constructed the objective function
M (x) = F (x, t) =

n
X
i=1

1
|targeti |

Z

t

(τ + 1)|ei (x, τ )| dτ,

(2.7)

t0

where n, targeti and 0 ≤ t0 < t are given and the error functions ei are
computed by the simulation software.
2.4. Character of the problem and practical limitations
Character of the problem plays a crucial role when choosing a fitting
method to solve it. For any set of controllers’ gains, the simulation provides
us with all data necessary to compute the objective function value (equation (2.7)). However, we cannot analyze the model itself and that forces us
to take it as a black box. We can only make a few careful assumptions.
In general, objective function M is non-convex, non-linear, non-quadratic
and highly multimodal (i.e. having multiple local optima). It is probably
continuous, but we know nothing of its conditioning (i.e. when a function
is ill-conditioned, it may be very steep in some places), and there are no
derivatives available. We cannot presume that M is smooth, so neither the
derivative approximations would necessarily be sensible. However, we do
not expect M to be noisy due to its heart made of a FEM model. We also
know that the objective function is non-separable, i.e. we cannot tune one
parameter independently of the others and expect to find the optimum that
way. One controller’s parameter’s value strongly affects other parameters
8

and one controller’s setting may affect the behavior of other controllers (the
controllers may be coupled). Since the relationships between the variables
are far from random, we would like the algorithm to mine and use this information to improve its search.
Beside this, we have to consider other practical limitations. Engine simulations, as in WAVE, take up to several minutes (per case – when parameters
like rotations per minute etc. are fixed), but the time budget is limited. Usually aiming for an overnight computation, we are allowed only a few thousand
evaluations of the objective function. That is, when parallelization is used,
and therefore we require the optimization algorithm to be parallelizable. Last
but not least, the algorithm must be robust – it must reliably produce decent
results, even at the cost of somewhat slower performance.
3. The optimization method
Metaheuristic and evolutionary methods have been extremely successful when tackling hard black-box optimization problems [6]. In order to
choose a method appropriate for our problem, we consider the extensive
Black-Box Optimization Benchmarking (BBOB) results, e.g. [20, 3, 4, 23],
which compare many algorithms on a wide range of test functions. Following the relevant data – tests on multimodal, non-separable, noiseless functions – we choose the Covariance Matrix Adaptation Evolution Strategy with
bi-population restart scheme (BIPOP-CMA-ES). Despite its efficiency and
growing number of practical applications, its use for tuning of PID controllers
is scarce [42, 25, 26].
3.1. CMA-ES: idea of the basic algorithm
The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [17,
16] is an evolutionary algorithm that uses stochastic and algebraic tools to
define optimally diverse population of candidate solutions in an area that
seems to be most promising. The size of the area and its location are determined based on the algorithm’s previous experience with the objective
function. New candidate solutions are sampled from a multivariate normal
distribution:
xk ∼ m + σN (0, C), k = 1, . . . , λ,
(3.1)
whose mean m and covariance matrix C are adapted in each generation along
with the general step size σ. The number of sampled candidate solutions λ
is called the population size.
9

The mean in the new generation is defined as weighted average of several best-ranking individuals (the parent set) of the last generation. The
weights are constant and depend only on the individuals’ relative ranking.
The covariance matrix is supposed to capture information about the objective
function’s features (curvature) around the mean. Detailed description of the
sophisticated updating of the covariance matrix can be found in the original
paper [17] or the informal tutorial [15]. In short, the matrix of the previous generation is weighted and further improved by exploiting the maximum
likelihood principle (when estimating the distribution’s parameters from the
data) and by adding information of the overall progress across generations
(the evolution path). The latter technique is also used for controlling the
overall step length σ.
Authors of CMA-ES have designed the method so that minimal user
interaction or manual setting of parameters is necessary. The starting point is
needed in order to initiate the first mean m(0) but otherwise all the parameters
are assigned usable default values. A user might want to change only the
initial step length σ (0) and the population size λ. The parent set size is given
as µ = 12 λ. The population size is by default very small, given by the formula
λdef = 4 + b3 log(d)c , where d is the problem dimension [17]. For a given
budget, smaller population size means faster adaptation. More often than
not, this is a desirable property. Choice of m(0) and σ (0) is described further.
3.2. Extensions of the basic algorithm
The basic algorithm described above can be further extended and upgraded. In our application, we shall use the elitist BIPOP-aCMA-ES version, i.e. Covariance Matrix Adaptation Evolution Strategy [17] with active
covariance matrix updates [27], elitist scheme of parent selection [21] and bipopulation restart strategy [22]. The choice of this particular variant, as well
as setting the algorithm’s parameters, is supported by numerical experiments
described in section 4.
3.2.1. Restart strategies
Restart strategies enable the basic algorithm to become more robust by
outweighing premature convergence. They are based on changing the population size, which also helps to balance exploration (exploring as much of the
vector space as possible) and exploitation (thorough investigation of promising areas) better.

10

The default population size λdef is small and the optimal population size
may be much larger. The IPOP (increasing population) strategy doubles
the population size (the factor of 2 is empirical) every time a new restart is
launched, while other parameters remain unchanged [2].
The more advanced BIPOP (bi-population) restart strategy [22] makes
use of two interlacing regimes. The first one uses the IPOP restart strategy,
while the second uses varying small populations. After each restart, it is
decided which of the two regimes is to be applied next depending on whose
count of conducted function evaluations is lower. The maximal population
size is limited by the number of restarts under the first regime.
There are two prominent restart criteria that need to be set carefully and
adjusted for each model, because they depend on its numerical values. The
first one is TolFunHist, which gives tolerance in function value history (range
of the best objective function values of several last generation is almost zero).
The second criterion is TolFun, which gives tolerance in function values of
the current generation (stop, if the difference between the current best and
worst candidate solution is TolFun-negligible) and it also checks the history
of function values in a very similar way to TolFunHist. Their default values
are too small for our use and better values can be guessed from the objective
function’s numerical values.
It was observed that it is generally good to try larger values of TolFunHist
first and diminish it if it invokes restarts too often. When logging the algorithm’s run, occurrence of premature restarts is usually obvious. Criterion
TolFun can save some runtime too, when it is set up correctly (in the same
way as the first criterion). Or we can set it to a small value, preferably a
fraction of TolFunHist. That way, most restarts are called by TolFunHist
and we need to tune only this one restart criterion.
3.2.2. Elitist selection
There are multiple options of how to choose parents of a new generation
[21]. The basic algorithm selects all parents from the current generation
(non-elitist selection). However, this scheme does not exploit very good solutions found early in the search. They are used as parents within the one
generation that produced them and then are lost. Selecting the parent set
among the individuals of the current generation and their parents as well
(elitist selection) helps to preserve the exceptionally good individuals until
they are superseded and thus amplify their influence. This approach speeds
up the convergence and is advantageous in our problem, where good solu11

tions are scarce. The downside is that it may, and sometimes it does, lead
to premature convergence to a local optimum.
3.2.3. Active covariance matrix adaptation
In the original algorithm, the successful individuals are used for the covariance matrix adaptation. Variance in directions that have proven to be
beneficial is increased and thus they are preferred when sampling next generation. With the so-called active updates, information hidden in the unsuccessful individuals is exploited as well [27]. As opposed to passive decay over
time, variance in detrimental directions is actively decreased. Simply put,
besides telling the method where to go, we also tell it where not to go.
3.3. Important properties of CMA-ES and their use
CMA-ES does not use gradients and it does not even presume their existence. Moreover, it does not even use the actual values of the objective
function once relative ranking has been assigned to the candidate solutions
(except for some stopping criteria). As a result, it is invariant to strictly
monotonic transformations of the objective function. I.e. those transformations of the objective function that have no effect upon the relative ranking of
individuals do not effect the method’s performance, making it more robust.
Further, the method exhibits invariance to invertible linear transformations of the search space. In particular, CMA-ES is invariant to scaling of
variables (coordinate axes), which is the key property that makes it wellsuited for tuning multiple controllers: parameters of one controller are usually of roughly the same scale, but with multiple controllers, the scaling
may differ by many orders. Scaling of variables is used in the following
way. The algorithm is given an initial approximation of the solution (hereafter called the reference point): s = (s1 , . . . , sN ). Then, when the algorithm wants to evaluate vector v = (v1 , . . . , vN ), it scales it by the reference point’s elements’ magnitudes and inputs vector w in the engine model:
w = (|s1 | · v1 , |s2 | · v2 , . . . , |sN | · vN ). Because the scale coefficients are all
positive and signs of the variables are essential, the algorithm must always
have a zero vector as its starting point. Also, the initial step length is set to

12

σ (0) = 1.0.
Algorithm 1: Elitist BIPOP-aCMA-ES
set λ, µ
initialize m, σ, C = I, pσ = 0, pc = 0
initialize restart regime = 1, count1 = 0, count2 = 0
while termination criteria not met do
while restart criteria not met do
if not first generation in a restart then
for i = 1, . . . , µ do
xi+µ = xi
// relabel parents of previous generation
fi+µ = fi
// relabel parents’ objective function values

for i = 1, . . . , λ do
xi ∼ N (m, σ 2 C)
fi = evaluate(xi )

// sample new population from normal distribution
// evaluate xi with objective function

sort xi , i = 1, . . . , λ + µ acc. to fi
// assign relative (descending) ranking
m∗ = m
m = update m(xi , . . . , xµ )
// move the mean utilizing the parents
// the evolution paths contain information about past progress
pσ = update pσ (pσ , σ −1 C −1/2 (m − m∗ ))
// isotropic evolution path update
pc = update pc (pc , σ −1 (m − m∗ ), kpσ k)
// anisotropic evolution path update
C = update C(C, pc , (x1 − m∗ )/σ, . . . , (xλ+µ − m∗ )/σ) // covariance matrix update
σ = update σ(σ, kpσ k)
// step size update
if restart regime = 1 then
count1 = count1 + λ
else
count2 = count2 + λ

if count1 < count2 then
restart regime = 1
else
restart regime = 2
reinitialize parameters and variables acc. to selected restart regime

4. Experiment-based calibration of the optimization method
For the experiments, the 1D engine simulation software package WAVE
by Ricardo Software was used [36], as the primary incentive was to develop
a working automated tuner of PID controllers within WAVE. For CMA-ES
implementation, we take the Python code distributed by Hansen [13].
Since models of engines take several minutes to run, while the CMA-ES
13

Figure 2: The basic testing model in WAVE.

runtime is negligible, we shall measure time in number of evaluations of the
objective function (= number of simulation runs).
4.1. The basic testing model
For the algorithms’ calibration and basic testing, a quick-to-run model
was used. It has three strongly coupled controllers and, despite its relatively
simple mechanics, it proved hard to tune.
In this model, we can see a single cylinder (orange circle) engine. The blue
“clouds” contain information about the surroundings (e.g. ambient pressure
and temperature or initial fluid composition). The thick black lines depict
the ducts. The green element is an orifice – an opening of variable diameter.
The yellow PID elements are the controllers, whose gains we want to tune.
From left to right, the first PID controls the orifice diameter, the second
one controls the fuel-air ratio by manipulating the fuel injector, and the
third controls the compression ratio (the ratio of the maximum to minimum
volume in the cylinder). The arrow-like elements are actuators that perform
the actual mechanical control based on the control signal outputted by the
corresponding controllers.
14

Figure 3: The basic testing WAVE model: comparison of good solutions. Pink: the target
value, dark blue: tuned by engineer, light blue: tuned by the algorithm.

The controlled quantities measured by sensors (depicted as gray circles)
are: indicated mean effective pressure (IMEP; the average pressure acting
upon the piston during its cycle; controlled by adjusting the fuel-air ratio),
exhaust gas temperature (controlled by the compression ratio, i.e. the ratio
of largest and smallest possible capacity of the combustion chamber), volumetric efficiency (the ratio of the volume of fluid actually displaced by a
piston; controlled by opening of the orifice). Clearly, they have nontrivial
influence on each other, so the controllers are coupled and cannot be tuned
independently of each other.
We then simulate the start of this engine, when the controllers face a
simultaneous step change of all three target values. The target values of the
controlled quantities are given to us (appropriately to the context). In figure
3, each plot contains two step responses. Controllers tuned by an engineer
(hereafter called the baseline solution) yield the dark blue curves that describe the development of the controlled quantities’ values. This solution is
sufficient but we can see that a solution found by our method (light blue)
is even better. More importantly, it can be reached with very little user
interaction, saving a lot of expert manual work.
However, most candidate solutions found in the tuning process are not
acceptable: oscillatory, convergent to a different value in one or more criteria,
or even divergent (see figure 4). In either case, it is important to observe the
development of the quantities’ values in order to determine the shift t0 used in
the objective function: we choose a time point just before the initially erratic
curve starts to follow a trend. For this model of overall simulation time 6
15

Figure 4: The basic testing WAVE model: unsatisfactory solutions still provide enough
information to estimate a fitting shift.

seconds, the best time shift found is t0 = 1.2 second. It seems fitting to set
the total simulation time such that the shift makes up the first 20%. That is,
when all the controlled quantities react with approximately the same speed
as in this case. In general, we need to compromise between the slowest and
quickest responses: cut out the uninteresting information but leave enough
to guide the search. Larger t0 means greater risk of leaving out important
information but possibly faster search.
Setting t0 > 0 is specific to our problem and certainly cannot be used in all
applications of PID controllers. Whenever large overshooting is an issue, e.g.
due to circuit breakers, this approach is not applicable. In our case, however,
we have an empirical guarantee that leaving out most of the transient part
(i.e. when 0 < τ < t0 ) does not lead to unacceptable solutions. For the
purposes of optimization, it is even desirable to ignore it as that reduces the
“noise” in the objective function caused by this information. The problem
then becomes easier to solve for the optimization method.
4.2. Tuning CMA-ES
The basic testing model was used to tune CMA-ES and other attributes.
For that purpose, we consider all three controllers to be only PI controllers
(i.e. with derivative gains set to zero) as is common in our application. Because the method uses randomness and each run is different, each test is run
ten times, with no change to the given reference point. Then the minimum,
maximum and average run times are computed. The same reference point
(i.e. scaling) of rather ‘poor quality” is used for all tests. When we compare
it to the vector of a good solution, the elements are off by two orders of
magnitude on average.

16

A run is terminated when the target objective function value is hit or the
budget of 12000 function evaluations is depleted. In each run, a different solution might be found, but setting the target value to 0.5 ensures high quality
of any of them (all the corresponding step responses are almost identical).
In practice, however, the algorithm will have to be stopped manually by the
user because the desired objective function value varies greatly from model
to model and the target value is difficult to estimate.
The results are summarized in table 1. The best setting found is always
compared to the setting, where one or two specified attributes were changed.
setting
CMA-ES variants

min

max

average

active, not elitist

644

10877

4145

elitist, not active

327

5361

2184

not active, not elitist

812

> 12000*

> 5334

adaptation method

TPA

926

> 12000*

> 4724

λ0 , µ0

λ0 = 2λdef , µ0 = λdef

1533

8812

3273

λ0 = 2λdef , µ0 = 1/2λdef

697

10973

4181

0.8 × 1.5

987

8980

3147

1.2 × 1.5

674

4990

2058

2.0 × 1.5

871

4960

2837

0s

919

9567

3687

0.6 s (10%)

354

11169

2963

268

2267

1098

TolFunHist

shift t0

BEST setting

Table 1: Experiments’ results comparing various settings of CMA-ES and the objective
function, where always only the specified attribute or two were changed from the BEST
setting. BEST setting: shift t0 = 1.2s (20%), CMA-ES variants = active and elitist,
adaptation method = CSA, population size λ0 = λdef , parent set size µ0 = 1/2λdef ,
restart parameter TolFunHist = 1.5. Entries marked *: number of evaluations once
exceeded the maximum number of iterations allowed. When computing the corresponding
average value, the value of 12000 was used in such cases, resulting in lower-bound for the
average run time.

When choosing a variant of CMA-ES, the results show that the best
combination is the elitist parent selection scheme with active updates of the
covariance matrix and CSA step size adaptation (the Cumulative Step-size
Adaptation [17] uses information of the algorithm’s overall progress across
generations). The alternative TPA (Two-Point step-size Adaptation, [14])
implements a line search along the direction of the latest mean shift.
17

Then we experiment with the initial population size λ0 (i.e. the population size before the first restart, which is further used to compute the
population sizes in the next ones). Using smaller than default populations is
discouraged by Hansen [15], while more exhaustive search with bigger populations motivated the development of the IPOP restart strategy. Therefore,
we try to double the initial population size: λ0 = 2λdef . The parent set size
µ0 , by default equal to half of the population size, is then either doubled as
well (i.e. µ0 = 12 λ0 ) or it stays the same (i.e. µ0 = 12 λdef ). For a given cost,
smaller populations enable more generations than large populations, causing
faster adaptation [17], which seems beneficial in our application.
The best working value of the influential restart parameter TolFunHist
(tolerance in function value history) for the basic testing model was found
to be n/2, where n = 3 is the number of controlled quantities. It was
further observed that setting n/2 seems to be almost universally usable in our
context, even though its optimal value differs for each model. The tolerance
in function value was set to TolFun = 0.1 for the basic testing model and
later roughly adjusted for other models, scaled in proportion to the “usual”
numerical values of the objective function and expected threshold.
In the objective function, shift t0 = 1.2 second proved to be most effective.
4.3. Testing robustness of the method
Now we compare the tuned method’s performance with various reference
points. We take the baseline solution and derive the other reference points
by multiplying each of its elements by factors 10−3 , 10−2 , 10−1 , 101 , and 102
(scaling by 103 proved to be too challenging for all tested methods). These
settings simulate user estimates of various quality.
In these tests, both PI and PID controllers are considered. Since we do
not have a PID baseline solution with nonzero D parameters, we estimate
them based on the corresponding magnitudes of P and I parameters.
Results are summarized in table 2 together with other methods’ robustness tests. We can see that the run times are very good: of 120 runs, 8
reached 2000 to 3000 evaluations and only two exceeded 3000 evaluations.
Surprisingly, there is no greater difference between the 6-dimensional PI setting and 9-dimensional PID setting. Higher problem dimension is balanced
by greater flexibility of the system. Closer look at the behavior of the algorithms suggests that too large search area (i.e. scaling) causes them to “get
lost” very far from the optimum. The run times are also influenced by the
fact that there exist good solutions at the “10−1 level”.
18

We can conclude that when the reference point (scaling) is within a reasonable range of two orders of magnitude in each coordinate, the average
run times are very usable. A user should be able to provide such a reference
point – an estimate of the solution. If in doubt, magnitudes of the reference
point’s coordinates should be chosen rather smaller than larger.
5. Comparison with PSO and SHADE
Performance of the above-described variant of CMA-ES on the problem
of tuning coupled controllers was further compared with performance of two
other prominent evolutionary algorithms: the particle swarm optimization
(PSO) and the differential evolution (DE) in its success-history based adaptive variant SHADE.
5.1. PSO
PSO [28, 34] is inspired by bird flocks or fish schools, where every individual moves by itself, yet the whole self-organized system acts as a single
organism. Each particle within a swarm moves in the search space as it is
assigned a different “velocity” vector in each generation. This vector is defined as a combination of the previous “velocity”, the individual’s best known
position so far and the swarm’s (or sub-swarm’s) best known position. This
way, the whole swarm moves towards historically best areas.
For PSO implementation, we use the default setting of PSO in the DEAP
(Distributed Evolutionary Algorithms in Python) optimization framework
[7]. We only add the essential scaling of variables same as in CMA-ES.
5.2. SHADE
The basic DE algorithm [35, 37, 38] maintains a population of candidate
solution, which are further combined and tested for better objective function
values. Each candidate solution X is tested against a new point Z, which
was obtained as a binary crossover of X and Y = A + f (B − C), where A, B
and C are three distinct candidate solutions in the population and f ∈ [0, 2]
is a parameter. The crossover probability of each vector element is given
by parameter cr ∈ [0, 1]. Both parameters and the population size greatly
influence efficiency of the search.
Considering the context of our problem and the performance-comparing
tests in [41], the SHADE variant of DE was chosen to compete with CMAES. The SHADE upgrade [39] introduces adaptation of the DE parameters
cr and f based on success history that is stored in an external archive.
19

20

56
893

102 PI b.

PI b.

10−1

10−2 PI b.

33
387

101 PID b.

102 PID b.

PID b.

10−1

10−2 PID b.
-

PID b.

calibration ref. p.

10−3

562

PID baseline

-

-

101 PI b.

PI b.

113

PI baseline

10−3

#1

reference p.

-

-

624

164

-

-

28

-

-

76

-

-

61

#2

-

1362

301

104

-

-

105

-

483

48

-

-

-

#3

280

-

-

71

-

226

39

-

-

-

#4

-

1734

1186

PSO runs

-

1680

-

43

-

-

28

1600

403

32

-

-

-

#5

-

1673

6524

239

-

-

428

4118

5155

342

-

-

130

#1

-

2125

1018

418

-

-

307

1090

1514

156

-

-

396

#2

-

1945

1150

209

-

-

497

2102

2783

1560

-

-

341

#3

-

1477

3297

437

-

-

441

1243

1437

260

-

-

2100

#4

SHADE runs

-

2163

1556

390

-

-

658

1422

2123

538

-

-

155

#5

268

416

202

141

953

62

32

259

222

61

821

540

12

min

2267

2324

941

782

4130

1462

256

1617

1334

816

5700

2064

168

max

CMA-ES runs

-

-

-

125

-

-

159

-

-

50

-

-

-

PSO

-

1877

2709

339

-

-

466

1995

2602

571

-

-

624

SHADE

1098

1138

580

343

2022

1102

67

812

592

317

2061

1049

76

CMA-ES

average run time

Table 2: Results of algorithm testing on the basic model. 5 PSO runs, 5 SHADE runs and 10 CMA-ES runs were performed
for each of 13 reference points. The value of “-” means that a satisfactory solution (i.e. solution with function value less than
0.5) was not found within the provided budget of 10000 function evaluations. Average run time was not computed if one or
more runs did not finish within the given budget.

Implementation of the SHADE algorithm by Tanabe was used [40] and
tuned in accordance with its author’s recommendations [41] and additional
experiments, setting, most notably, initial cr = 0.5, initial f = 0.5 and population size = 2N , where N is the problem dimension. One restart criterion
was added to the standard SHADE algorithm, based on recommendations
in [41]: if the number of objective function evaluations exceeds 1000, restart
can be launched if the best-so-far solution is not updated for 50 × N evaluations. The particular numbers were set based on experiments and the first
condition proved very important.
Unlike in CMA-ES and PSO, the reference point is not used for scaling
but to define boundaries of the area, from which the initial population is
sampled. Each coordinate of the reference vector is multiplied by −10 and
10 to define the boundaries. Candidate solutions in further generations are
not restricted.
5.3. PSO and SHADE experiment results
The experiments with PSO and SHADE were performed with same set of
13 reference points as in the previous tests of CMA-ES, with 5 runs of each
test. The stopping criteria were either hitting the target value of the objective
function or depleting the budget of 10000 objective function evaluations. All
test results can be found in table 2.
PSO was able to outperform CMA-ES but was not very reliable. Its good
performance was limited to very good reference points and the algorithm very
often (36 out of total 65 runs) did not converge within the given budget.
The performance of SHADE was more consistent but very slow compared
to CMA-ES. Same as PSO, SHADE never converged within the given budget
for the reference points of baselines multiplied by 101 and 102 or the one used
for CMA-ES calibration.
6. Verification of the method on models of real engines
We have shown that, on the testing model, the described version of CMAES performs very well and clearly outperforms PSO and DE. Now we verify
its usability on models of real engines provided by Ricardo. The testing set
consists of two models with 1 controller (labeled M1.1, M1.2), three models with 2 controllers (labeled M2.1, M2.2, M2.3) and one model with 3
controllers (labeled M3.1). For each, we are given a baseline setting of controllers: a setting tuned manually by an engineer, which is always in the PI
21

model

M1.2

min

max

aver.

28

PI baseline

11

66

35

101

244

280

255

4

32

21

min

max

aver.

PI baseline

2

68

101
M1.1

reference p.

reference p.

PI b.

35

153

79

102 PI b.

95

519

225

10−1 PI b.

20

120

66

10−2 PI b.

49

296

123

PI baseline

1

22

9

101

PI b.

4

28

11

102 PI b.

80

225

187

10−1 PI b.

34

100

51

10−2 PI b.

57

181

94

model

M2.1

PI b.

10−1 PI b.

M2.2

PI baseline

8

98

29

101 PI b.

60

770

364

10−1 PI b.

44

107

64

PI baseline
M2.3

M2.3

9

78

32

101 PI b.

250

757

629

10−1 PI b.

49

1188

347

PID baseline

10

91

57

101

274

857

522

82

1576

749

PID b.

10−1 PID b.

M3.1

PID baseline

41

331

152

101 PID b.

827

1763

1268

10−1 PID b.

179

3867

2476

Table 3: Real-world models

form (i.e. with D gain set to zero). When a baseline value of the D gain is
needed in our tests, we estimate its magnitude from the P and I gains.
The single-controller models are tested with 5 reference points: the baseline solution and its element-wise multiples by factors 10−2 , 10−1 , 101 and 102 .
The multi-controller models are tested with three reference points: the baseline solution and its multiples of 10−1 and 101 . The primary reason for this
restriction is the enormous time consumption. Each test is run five times. A
run is terminated upon reaching the (empirically set) target value specific to
the given model. Even though each run might produce a different solution,
as long as the target value is reached, their quality is guaranteed.
All results are summarized in table 3.
6.1. Single-controller models
Even though our primary objective is to tune multiple controllers, we
include also single-controller models as they are very important in practical
use and we want our method to work for them as well. Model M1.1 contains
a turbocharger and M1.2 contains a twin turbocharger. They represent the
typical use of a PI controller in a combustion engine.
22

With the model dimension being only 2, the run times for both M1.1
and M1.2 are very short. Good reference points for single controllers can be
obtained by commonly used rules of thumb (e.g. Ziegler-Nichols) or provided
by previous (personal or programmed) experience with similar models.
6.2. Multi-controller models
The two-controller models M2.1 and M2.2 do not pose a greater challenge
than the single-controller models. For model M2.3, we compare PI and PID
control. The extreme differences in minimal and maximal run time values
of “10−1 PI baseline” and “10−1 PID baseline” tests are caused by nonoptimal setting of the restart parameter tolerance in function value history
(tolhistfun) described above. After adjustment, the “10−1 PI baseline” and
“10−1 PID baseline” run times drop and the results of PI control resemble
those of M2.2.
Unlike the previous cases, the three-controller model M3.1 requires a full
PID control. The reference point was obtained by estimating the D gains
and adding them to the given (not very good) PI engineer-tuned baseline
“solution”. It can be seen that this model is considerably harder to tune
than the previous models. While the algorithm has no trouble finding the
same near-optimal “solutions” similar to the one given by an engineer, it was
hard to get to an actual global optimum, when all three controllers converge.
7. Conclusion
This paper has shown how the Covariance Matrix Adaptation Evolution
Strategy can be applied to the problem of tuning the gains of multiple coupled
PID controllers within combustion engine simulations. We have shown that
its version with bi-population restart scheme, elitist parent selection and
active covariance matrix updates is capable of finding good parameters of up
to three PID controllers through minimization of a fitting objective function.
The method has been calibrated on a testing model and verified on models of
real-world engines, showing its practical usability and tolerable computation
times even for poor-quality reference points. On the testing model, CMA-ES
clearly outperformed PSO and SHADE methods.
Acknowledgments
This research did not receive any specific grant from funding agencies in
the public, commercial, or not-for-profit sectors. The author would like to
23

thank Ricardo Prague, s.r.o. for their support and to anonymous reviewers,
whose comments and suggestions inspired significant improvements.
References
References
[1] K. J. Åström and B. Wittenmark. Adaptive Control: Second Edition.
Courier Corporation, 2013. ISBN 9780486319148.
[2] A. Auger and N. Hansen. A restart CMA evolution strategy with increasing population size. In B. McKay et al., editors, The 2005 IEEE
International Congress on Evolutionary Computation (CEC’05), vol. 2,
pp. 1769–1776, 2005.
[3] A. Auger, S. Finck, N. Hansen, and R. Ros. BBOB 2009: Comparison
Tables of All Algorithms on All Noiseless Functions. Technical Report
RT-0383, INRIA, April 2010.
[4] A. Auger, S. Finck, N. Hansen, and R. Ros. BBOB 2010: Comparison
Tables of All Algorithms on All Noiseless Functions. Technical Report
RT-388, INRIA, September 2010.
[5] Z. Bingul. A new pid tuning technique using differential evolution for unstable and integrating processes with time delay. In Neural Information
Processing: 11th International Conference, ICONIP 2004, Calcutta, India, November 22-25, pages 254–260, Berlin, Heidelberg, 2004. Springer
Berlin Heidelberg. ISBN 978-3-540-30499-9.
[6] J. Dréo, A. Petrowski, P. Siarry, and E. Taillard. Metaheuristics for hard
optimization: methods and case studies. Springer Science & Business
Media, 2006.
[7] F.-A. Fortin, F.-M. De Rainville, M.-A. Gardner, M. Parizeau, and C.
Gagné. DEAP: Evolutionary Algorithms Made Easy. Journal of Machine Learning Research, 13:2171–2175, jul 2012.
[8] Z.-L. Gaing. A particle swarm optimization approach for optimum design of PID controller in AVR system. IEEE Transactions on Energy
Conversion, 19(2):384–391, June 2004. ISSN 0885-8969.
24

[9] S. P. Ghoshal. Optimizations of PID gains by particle swarm optimizations in fuzzy based automatic generation control. Electric Power
Systems Research, 72(3):203–212, 2004.
[10] D.E. Goldberg. Genetic Algorithms in Search, Optimization, and Machine Learning. Artificial Intelligence. Addison-Wesley Publishing Company, 1989. ISBN 9780201157673.
[11] G. C. Goodwin, S. F. Graebe, and M. E. Salgado. Control System
Design. Prentice Hall PTR, Upper Saddle River, NJ, USA, 1st edition,
2000. ISBN 0139586539.
[12] D. Graham and R. C. Lathrop. The Synthesis of Optimum Transient
Response: Criteria and Standard Forms. Transactions of the American
Institute of Electrical Engineers, Part II: Applications and Industry, 72
(5):273–288, Nov 1953. ISSN 0097-2185.
[13] N. Hansen. cma 1.1.06: Python Package. https://pypi.python.org/
pypi/cma.
[14] N. Hansen. CMA-ES with Two-Point Step-Size Adaptation. Research
Report RR-6527, INRIA, 2008.
[15] N. Hansen. The CMA Evolution Strategy: A Tutorial. https://www.
lri.fr/~hansen/cmatutorial.pdf, 2015.
[16] N. Hansen and S. Kern. Evaluating the CMA evolution strategy on multimodal test functions. In International Conference on Parallel Problem
Solving from Nature, pages 282–291. Springer, 2004.
[17] N. Hansen and A. Ostermeier. Completely derandomized self-adaptation
in evolution strategies. Evolutionary Computation, 9(2):159–195, 2001.
[18] N. Hansen, A. S. P. Niederberger, L. Guzzella, and P. Koumoutsakos.
Evolutionary optimization of feedback controllers for thermoacoustic instabilities. In IUTAM Symposium on Flow Control and MEMS, pages
311–317. Springer, 2008.
[19] N. Hansen, A. S. P. Niederberger, L. Guzzella, and P. Koumoutsakos. A
method for handling uncertainty in evolutionary optimization with an
application to feedback control of combustion. IEEE Transactions on
Evolutionary Computation, 13(1):180–197, 2009.
25

[20] N. Hansen, A. Auger, R. Ros, S. Finck, and P. Posik. Comparing Results of 31 Algorithms from the Black-Box Optimization Benchmarking
BBOB-2009. Workshop Proceedings of the GECCO Genetic and Evolutionary Computation Conference 2010, pages 1689–1696, 2010.
[21] N. Hansen, D.V. Arnold, and A. Auger. Evolution Strategies. In J.
Kacprzyk and W. Pedrycz, editors, Springer Handbook of Computational Intelligence, chapter 44, pages 871–898. Springer Berlin Heidelberg, 2015. ISBN 9783662435052.
[22] N. Hansen. Benchmarking a BI-population CMA-ES on the BBOB2009 Function Testbed. In Proceedings of the 11th Annual Conference
Companion on Genetic and Evolutionary Computation Conference: Late
Breaking Papers, GECCO ’09, pages 2389–2396, 2009. ISBN 978-160558-505-5.
[23] N. Hansen and S. Finck and R. Ros. COCO - COmparing Continuous
Optimizers : The Documentation. Research Report RT-0409, INRIA,
May 2011.
[24] J. H. Holland. Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial
intelligence. University of Michigan Press, 1975.
[25] M. W. Iruthayarajan and S. Baskar. Evolutionary Algorithms Based
Design of Multivariable PID Controller. Expert Syst. Appl., 36(5):9159–
9167, July 2009. ISSN 0957-4174.
[26] W. M. Iruthayarajan and S. Baskar. Covariance Matrix Adaptation
Evolution Strategy Based Design of Centralized PID Controller. Expert
Syst. Appl., 37(8):5775–5781, August 2010. ISSN 0957-4174.
[27] G. A. Jastrebski and D. V. Arnold. Improving evolution strategies
through active covariance matrix adaptation. In IEEE Congress on
Evolutionary Computation – CEC 2006, pages 2814–2821, 2006.
[28] J. Kennedy, R. C. Eberhart and Y. Shi. Swarm Intelligence. Morgan
Kaufmann Publishers Inc., San Francisco, CA, USA, 2001. ISBN 155860-595-9.

26

[29] W. M. Korani, H. T. Dorrah, and H. M. Emara. Bacterial foraging
oriented by Particle Swarm Optimization strategy for PID tuning. In
Computational Intelligence in Robotics and Automation (CIRA), 2009
IEEE International Symposium on, pages 445–450, Dec 2009.
[30] D.P. Kwok and F. Sheng. Genetic algorithm and simulated annealing
for optimal robot arm PID control. In Evolutionary Computation, 1994.
IEEE World Congress on Computational Intelligence., Proceedings of
the First IEEE Conference on, pages 707–713. IEEE, 1994.
[31] R. T. Marler and J. S. Arora. Survey of multi-objective optimization
methods for engineering. Structural and multidisciplinary optimization,
26(6):369–395, 2004.
[32] K. Ogata. Modern Control Engineering. Prentice Hall PTR, Upper
Saddle River, NJ, USA, 2nd edition, 1990. ISBN 0135891280.
[33] M. Petrov, I. Ganchev and A. Taneva. Fuzzy PID control of nonlinear
plants. Intelligent Systems, 2002. Proceedings. 2002 First International
IEEE Symposium, vol. 1, pp. 30–35. IEEE, 2002.
[34] R. Poli. Analysis of the publications on the applications of particle
swarm optimisation. Journal of Artificial Evolution and Applications,
2008:3, 2008.
[35] K. Price, R. M. Storn and J. A. Lampinen. Differential Evolution: A
Practical Approach to Global Optimization. Springer-Verlag New York,
Inc., Secaucus, NJ, USA, 2005. ISBN 3540209506.
[36] Ricardo Software. WAVE Manual, 2016.
[37] R. Storn. On the usage of differential evolution for function optimization.
In Fuzzy Information Processing Society, 1996. NAFIPS., 1996 Biennial
Conference of the North American, pages 519–523, Jun 1996.
[38] R. Storn and K. Price. Differential evolution – a simple and efficient
adaptive scheme for global optimization over continuous spaces. Technical Report TR-95–012, International Computer Science Institute, Berkeley, 1995.

27

[39] R. Tanabe and A. Fukunaga. Success-History Based Parameter Adaptation for Differential Evolution. IEEE CEC, 2013, pp. 71–78.
[40] R. Tanabe and A. Fukunaga. Evaluating the performance of SHADE
on CEC 2013 benchmark problems. IEEE CEC, 2013, pp. 1952–1959.
[41] R. Tanabe and A. Fukunaga. Tuning Differential Evolution for Cheap,
Medium, and Expensive Computational Budgets. IEEE CEC, 2015, pp.
2018-2025.
[42] Y. Wakasa, S. Kanagawa, K. Tanaka, and Y. Nishimura. Pid controller
tuning based on the covariance matrix adaptation evolution strategy.
IEEJ Transactions on Electronics, Information and Systems, 130(5):
737–742, 2010.
[43] J. G. Ziegler and N. B. Nichols. Optimum Setting for Automatic Controller. ASME transaction, pp. 759–768, 1942.

28

