Decision Engineering of an Decyl-lactose Generation Chemical
Laboratory Process Assisted by a Diploid Genetic Algorithm and a
Multicriterion Aggregation Method
S. MASSEBEUF', C. FONTEIX’, L. N. KISS”, I. MARC!
' Laboratoire des Sciences du Génie Chimique —- UPR CNRS 6811
ENSIC, | rue Grandville, BP 451, 54001 Nancy Cedex, France
> Université Laval — Faculté des Sciences de I’ Administration
Bureau 2537/PAP, Sainte-Foy G1K 7P4, Québec, Canada
Abstract. In many, if not most, optimization problems, industrialists are often confronted with
multiobjective de n problems. For example, in manufacturing processes, it may be necessary
to optimi: veral criteria to take into account all the market constraints. So, the purpose
choose the best tradeoffs among all the defined and conflicting objectives. In multicriteria
zation, after the decision maker has chosen all his objectives, he has to determine the
feria optimal zone by using the concept of domination criterion called the Pareto
domination. Two points, in the research domain, are compared. If one is better for all attributes, it
is a non dominating solution. All the non dominating points form the Pareto’s region. In this
paper, a multiobjective optimization algorithm is used to obtain this zone, based on a diploid
genetic algorithm. This is compared to two industrial applications: food’s granulation and
decyllactose synthesis. In the optimal zone, the decision maker has to choose the best solution
after he has made a ranking with all potential solutions. Finally, a Decision Support System shell
is developed in order to classify all solutions.
1 INTRODUCTION
In a lot of domains, processing or product formulation depend on several objectives to
take into account their different features. Extrusion processes for food are a good illustration for
this case. For example, it may be necessary to optimize different parameters such as texture,
flavour, and so on, in order to formulate a new product and to maximize yields, production and
product quality with minimal investment whereas these criteria are not optimal for the same
working conditions. In the engineering processes domain, multiple objectives have usually been
combined, often through a linear combination, to form a scalar objective function [8] or a
geometric average of these ones [17]. Another classical method consists in optimizing one
chosen answer and in turning into constraints the others [14]. These techniques depend on the
user’s choice at the begining, so preferences can bias the results.
is
Methods incorporating a domination criterion are often more interesting because they are
more general, more accurate and without a priori knowledge. The object is to find a non
dominated zone in which a decision maker will be able to choose the best solution. This region,
called the Pareto domain, is the set of all non dominated points. For the Pareto’s domination
criterion, a point dominats another if it is better for one criterion, and better or equal for the
others. The aim of the study is not to find immediately the preferred solution but to exclude all
conditions which are not interesting. This information allows to choose the industrial decision
maker’s preferences, because it restricts all possible choices, mistakes and bias.
Schaffer [15] describes an extension of the traditional Genetic Algorithm (GA) which
allows the searching of the parameter space where multiple objectives are to be optimized. His
VEGA (Vector Evaluated GA) gives a selection preference to the non dominated members of a
population. But, only extreme points on the Pareto front seem to be found with VEGA, so Horn
et al. [9] apply a niching pressure to spread its population out along the Pareto optimal tradeoff
surface. Fonseca and Fleming [12], in the same way, use the non domination ranking for
selection to move a population towards the Pareto front and a niching mechanism, such as
1
sharing [13], to keep the GA from converging to a single point. Moreover, a direct intervention
of an external Decision Maker (DM) gives interactive information in the multiobjective
optimization loop. So, a satisfactory solution of the problem is found as soon as the knowledge
is acquired [12].
In this study, the Pareto domain is obtained at first without external influences. The
algorithm use an adapted GA, and is applied to several industrial. On the last step, the DM
influences are taken into account and notions of decision engineering are required. So, the global
multicriterion analysis leads to the best tradeoff.
2 THEORITICAL ALGORITHMIC BACKGROUND TO THE MULTIOBJECTIVE
OPTIMIZATION
2.1 An improved Diploid Genetic Algorithm (DGA)
The multiobjective optimization methods developed use a DGA whose principles were
elaborated by Perrin ef al. [4].The diploid version is kept because its performances were found
to be better compared with a haploid one [6]. Each individual (which can be a possible solution
of the problem) is described by a four-tuple (aj, a;’, Dj, xj). a; and a;’ represent the two alleles of
the j gene, Dj is the dominance of one allele over the other chosen in {0,1} values and x;
represents the phenotype which is the result of the combination of the respective alleles, aj and
aj’:
xj = Dj.aj + (1-D)).a°
An initial population is created by generating a set of m points from the search area. Each
point is tested and evaluated. If this population is not the solution, then genetic operators are
used to make it evolve. Only the better individuals will survive (elitist selection) and participate
in the creation of a new generation. The reproduction of the individuals in the diploid model
consists of a multi-crossover of the two chromosomes of each parent, a mutation and a
homozygosity. Mutation consists, for the selectionned individual, of a randomly draw for all the
genes of the two chromosomes and all the dominances. Homozygosity allows to modify a child
by copying the phenotype x of the two chromosomes a and a’ (x; = aj = aj’). All the genetic
operators and their usefulness are detailed in [4]. If the generated child is worse than the worst
parent, he is not adapted. He is eliminated and another is created. So dominance may be
controlled to decrease the child death rate.
The dominance may be randomly chosen in {0,1} values. But, to keep dominance
features, table 1 shows how to control Dj:
Dominance of aj in the first parent
Recessive Dominant
Dominance of a;’ in Recessive Random D=l
The second parent Dominant D,=0 Function of the best parent
Table 1: choice of the value of the dominance Djé {0,1}
Statistically, children have more chances of being better than the parents with the control of the
dominance, so the convergence of the algorithm speeds up and less points are tested.
The population of each generation is evaluated until it satisfies the stop criterion:
fax — fmin < € where fmin and fmax are respectively the minimal and the maximal objective
function values in the current population and € is the given precision for solution estimation.
2.2 Use of one function to define the set Pareto
Previously, a multicriteria optimization algorithm was elaborated by Viennet ef al. [7],
but it need three steps, and the number of points in the Pareto set was not controlled. The aim of
this section is to define one function to make a monocriteria optimization with the DGA. A set
of m given points is randomly generated. All individuals are evaluated by the calculation of each
objective fj (i = 1,...n). Then, for each point x;, a value is associated, F(x;), the number of times
that it is dominated by all the others in the current generation:
F(x,) = > By where hip = | if xp dominate x;
p=!
hip = 0 else
Let be denoted by n the number of criteria. For two points x; and Xp in the same population, and
for all criteria i (i = 1,..n):
if f\(x)) WORSETHAN fi(Xp) Cijp = 0
if f\(x)) BETTERTHAN fi(Xp) Cijp =n + 1
if f\(x)) EQUALTO fi(xp) Ciip = 1
where cijp is an intermediate variable.
In the case of minimization WORSE THAN is equivalent to <, and BETTER THAN to >.
if Vice <n then hip = | and hy; = 0 because x; is dominated by xp.
isl
if Yc, Zn and Yc, 2n then hip =hy=0 (hj =0).
ist ist
So, for each point, a value F corresponds by applying the Pareto’s domination criterion.
Then, all the points may be classified. The value of the function for the better individuals in the
current generation is equal to zero. The classified individuals and their function value are
presented in Table 2 in the k”* generation:
xX xX) x2 Xe erence Xs Xex comme Xt Ken amene Xa
F 0 0 0 sane 0 F(R) cassis FOG) F061) cen F(X)
Table 2: classified individuals in the k" generation
The s better individuals are non dominated in the current population (k) and are the chosen
parents for the next (k+1). The (m-s) worse individuals are eliminated for the next generation
but are used for the children evaluation. Each new individual is compared with all m points in
the previous population. Then, a threshold t is defined to keep the better children:
t=s+E[0.3*(m-s)] where E[x] is the entire part of x
and 0.3 is empirically choosen
So, a child death rate is applied to the convergence when the new individual x; is worse than the
X, point, i.e. F(x;) > F(x,). When the new population is reformed, all the m points are evaluated
with the new values of the function F. The population of each generation is evaluated until the
stop criterion is satisfied: all the points are non dominated, that is to say that F is equal to zero
for the m points.
The DGA is used only once and the algorithm converges rather rapidly (3 to 5
generations, depending on the problem). The most interesting thing is that the number of points
in the Pareto’s zone is controlled whatever the sproblem studied. If 5000 points are required to
clearly draw the Pareto’s area, it is represented by all these points with a lot of accuracy and
without an increasing calculus time.
3 PARETO DOMAIN DETERMINATIONS
3.1 Extrusion process
Food granulation for cattle is presented as an industrial application. The goal is the
optimization of the working conditions of an extrusion process described by Courcoux et al. [8].
A pulverulent product is converted into granules due to the conjugated effects of heat, moisture
and pressure. The industrialist would like to simultaneously minimize moisture and friability of
his product and the energetic consumption of the process. For modelling these criteria, an
experimental design is realised [8]. Factors having an influence on this process are the feeding
rate, speed of rotation for mixing, flour temperature, speed of rotation and drawplate profile. For
this study, only two factors are taken into account, the flour temperature T, situated between 35
and 75 °C, and the drawplate profile D, between 2 and 6 cm (figure 2).
2
ao, 6
Ty Mt
Figure 2: The evolutions of energetic consumption, friability index and moisture of granules vs
the temperature and drawplate profile, by Courcoux et al. [8]
The optimal Pareto’s zone, obtained by our algoritm, is represented in figure 3 by black
points.
15
65
TPC) 55
45
yg Cosa tise tiiiitiiiy
2 3 4 5 6
D (cm)
Figure 3: representation of the optimal Pareto’s for the example of an extrusion process
The method give satisfactory results for the Pareto’s zone, and keeps 5000 points for its
representation. It allows to have a precise area with a controlled number of points. So, the
decision maker will be able to choose in this zone the best point.
3.2 Decylaction of lactose synthesis
The technique for Pareto region determination is applied to the definition of an optimal
working zone for the synthesis of a biosurfactant derived from lactose: the decyllactose. For this
4
study, 3 factors are considered, lactose concentration between 10 and 320 g/l, sodium
hydride/lactose molar ratio between 1 and 8 and decyliodide/lactose molar ratio between 1 and
4. Three responses are maximized: reaction yield (decyllactose produced/initial lactose
concentration), lactose conversion rate and productivity (decyllactose produced/time of
reaction). A fourth one is minimized: biosurfactant cost price. This reaction is carried out during
3 hours at 35°C.
In figure, triangles show positions of experiments obtained from Fuzzy Dynamic
Experiment Design, a technic developped in our laboratory. The banned and undefined zones
correspond to experimental conditions where product formation is impossible or improbable due
to physico-chemical conditions. Radial Basis Functions neural networks are used for modelling
each response.
molar ratio
a décyl/lactose
cm
2,54 } ih
2 al
QO
the
0 50 100 150 200 250 300 350
[lactose] g/l
molar ratio
egBLVlaCtO Se ts
7A
5 BAG
\0
4
a | oO
2
aT
0 50 100 150 200 250 300 350
[lactose] g/I
Figure 6 : Pareto set and experiments performed.
4 MULTICRITERION ANALYSIS FRAMEWORK
4.1 An original approach
Many real problems can be modeled as multicriterion problems but the decision
engeneering knowledge intervention in a specific chemical engineering competence field as the
cattle-food granulation process is a very interesting new research. We dare to believe that even if
our results are inscribed in a relatively narrow axis and in spite of their limits, they will
contribute to the development of the knowledge in the chemical process engineering and in the
decision making sciences.
This process control problem can obviously be most appropriately structured around an
ELECTRE type outranking synthesis approach (OSA) [11]. The main reasons justifying our
choice are, in fact, the realistic way in which this OSA links the preferences; this is somewhat
similar to the "democratic" principle much more natural when faced with a decision having
product quality and product profitability impacts.
4.2 Algorithmic design
Let us consider the input parameters’ set of the granulation process containing the
discretized temperature (between 35 and 75 C°) and drawplate profile D (between 2 and 6
cm) values on the one hand and the output parameters'set of the granulation process containing
three criteria, namely the energetic consumption, the friability index and the moisture of the
granules on the other hand. The criteria role has been attributed to the output parameters'set of
the granulation process, i.e. energetic consumption, friability index and moisture of granules (to
minimize the 3), even when the alternatives role has been attributed to the input parameters'set
of the granulation process, i.e. temperature and drawplate profile.
Our multicriterion analysis algorithm is essentially a comparison process where the
decision maker's preferences are expressed by means of some constraint-like parameters:
indifference, preference, veto thresholds and the weights of the criteria. In order to operationally
run the multicriterion analysis process (MAP), it is necessary to introduce the MAP's input set
‘Y with the following synoptic structure [3]:
B= TE, Lm, 1 Mya cajs Kjajs Pjay> Rays Qo} P fal-V fol W tat]
where
eT, - MAP's algorithm;
e m - number of alternatives;
e on - number of criteria;
© Moénxay ~ realm x n cardinal matrix containing the alternatives'evaluations in
relation to the criteria and m= f(xie Mimxal> L=L,...7myj = 1,...30j
e Kj,) - binary n cardinal vector indicating the nature of the criteria evaluation, with
| 0 if the evaluation in relation to the j" criterion is an ordinal type
kK,=
J
1 if the evaluation in relation to the j” criterion is a cardinal type
and k,e Kj,,) ; j= L...,n (in this paper all the criteria are cardinal type);
¢ ©) - binary n cardinal vector indicating the nature of the orientation of the criteria, with
b= 0 if the j" criterion is to be minimized
; 1 if the j® criterion is to be maximized
J
and 9, € ®j,);j = 1,....n (0 for all criteria in this paper);
e R,,) - real n cardinal vector containing the range for every criteria
and Pi€ Rap i eee
¢ Qj.) - teal n cardinal vector containing the indifference thresholds
and q ; (3 Qin d =Lesaht
e Pj,) - real n cardinal vector containing the preference thresholds
and pje Pasi 1 era
e Vj,) - real n cardinal vector containing the veto thresholds
and v , € Vind © Lyssstif
e W,,) - realn cardinal vector containing the criteria's relative importance coefficients w ;
and wie Wiad cm Weweys
The following relations must be verified between the thresholds:
0 $q,<p,<v,s Pyj=l.n
a
The relative importance coefficients are normalized, i.e. »y wj= 1.
j=l
This input set ‘ is processed using the MAP algorithm T, to determine the ranking from the
best to the worst (i.e. "nadir") alternative with possible ex aequo.
Cc [i il is calculated as follows:
cli: i] = y w,- fii]; ist,...m; =1,....m; it isj=d...n.
j=]
1 if-A,< q,
A, +
where 6, [.i]= D, 7 if q)<-A,Sp,
0 if p, <-A,
ym; i# 1. If some criteria call for
and A, = m,, - m;€ Minxajs J= 1
maximization while others require minimization, we have introduced two temporary sentries
(i.e. guards) variables k and k in our calculations of Aj.
If ¢, = 0, meaning that this criterion is to be minimized, then k = i and k =i, else if ¢, = 1
(criterion to be maximize) k = i and k = 1; this way A, [i]
1,...n;k=1,....m; k =1,...mj:k #k.
A discordance index D, [isi
level as with ELECTRE III:
0 if-A,< p
-_ A; om P;
[i:i]= Vi ~Pj
1 if v, <-A, t Shen tite hLagiti=lagM) Tl. ahi -i «
[ii]
for every pair of alternatives. These outranking degrees are obtained using the following
il = cliff) =
= M™.- Mp3 VOs5=
is also calculated for every pair of alternatives at each criterion
j
" if q,<-A,Sv,
Using the concordance and discordance indexes, we generate outranking degrees 0
camsisl,...m; i# i [10].
Rousseau-Martel formula: o
The obtained outranking relation sets may be represented as an outranking matrix.
Working with outranking relations implies synthetizing that matrix in such a way as to provide a
process control recommendation in the form of an alternative ranking.
We make that synthesis following the notion of outgoing flow and incoming flow from
PROMETHEE [16].
For calculations we use:
of= Yollsil ifi#i and o; = 2
1 if isi 1
and we get the total ranking of the alternatives from the net flow as in PROMETHEE II:
0, =0; -o; i=1,...,.m.
fsi)iviei_ ett tno it
1 if isi
4.3 Results
4.3.1 Extrusion
In order to illustrate the proposed treatment, firstly we present the results of the extrusion
process, using the earlier described I, algorithm. We have developed a Decision Support
System (DSS) shell based on our multicriterion analysis conception with the tresholds indicate
in table 3.
ch Pj Yj
Friability index 0.2 0.5 0.8
Moisture 0.5 1.5 3
Energetic consumption 1 3 6
Table 3: indifference, preference and veto tresholds for each criterion
Table 4 gives the best and the worst alternatives in the Pareto domain, after the total ranking has
been made:
Best Nadir
D (cm) TCC) D (cm) TCO
3.3186 64.909 5.2921 72.677
Table 4: best and worst parameters alternatives
This DSS shell produces very encouraging final results as process control recommandation for
the chemical engineering. To illustrate the robustness of the results, we present quintile by
quintile the input parameters of the granulation process as well as the positions of the best and
the worst input parameters (Figure 5). Each quintile contains 20% of the points in the Pareto set.
7 o— i -9
65 4 o |
TCC) 35 mY 4
ae
Figure 5: Best (#) and worst (m) conditions giving by DSS shell
We can notice that the best point is a robust one and the little variations of the working
conditions do not alter much the product quality, the point stay in the Pareto domain. On the
contrary, the worst point is on the borderline of the zone and so is not robust. The quintiles are
represented by concentric zones (first zone: », second zone: e, third zone: °, etc...).
4.3.2 Decylactation of lactose
Secondly we present the results concerning the decylactation of lactos process, using the
same Ty algorithm, but with the readjusted indifference, preference and veto tresholds for each
criterion applied in the decylactation of lactose context.
Table 5 gives the best and the worst alternatives in the Pareto domain, after the total ranking has
been made:
Best Nadir
Lactose (g/l) Molar ratio NaH/lactose lactose (g/l) Molar ratio NaH/lactose
97.889261 6.334321 115.52411 7.233606
Table 5: best and worst parameters alternatives for the decylactation of lactose process
The following figure illustrate the very encouraging robustness of the results, by presenting the
input parameters of the decylactation of lactose process as well as the positions of the best and
the worst input parameters (Figure 6).
8 +
Molar ratio NaHAactose
1 + + + + + +
a 50 100 160 200 760. 300 360
[lactose] (g.[")
Figure 6: Best () and worst (m) conditions giving by DSS shell
5 CONCLUDING REMARKS AND FURTHER CONSIDERATIONS
In this paper, a multiobjective decision problem is treated in order to explain an extrusion
process and an decyl-lactose generation chemical laboratory process The main objective is the
determination of the discussed chemical process's operating conditions which correspond not
only to the traditional optimization problems, but also to the choice of the operating conditions
according to the desired technological, chemical engineering priorities or preferences of the
decision maker. So, the proposed steps (at first, the determination of the multiobjective optimal
zone and then the multicriterion analysis) give a lot of data to choose the best tradeoff of the
problem.This type of analysis design concerns principally the decisional engineering
competence, but requires also a good knowledge of the chemical phenomena as well as the
definition of a sufficiently precise model.
In order to respect the paper length constraint, we precise that our summerizing remarks
regarde first of alls the granulation process, although one could formulate semantically
similar remarks also relating to the decylactation of lactose process.
Our DSS shell gives the complete ranking of the Pareto domain. We notice the
robustness of the best operating conditions, i.e. it can stay in the same quintile even with small
variations of temperature and profile drawplate. It is an important result for the decision making
science.
The results of our analysis can be registered in the chapter of the innovating contributions
of the art of chemical engineering and opens promising perspectives towards the exploration for
the decision robustness problem concerning the effective control of the chemical processes. For
example, one of the close perspectives is to apply the analysis to polymerization processes which
are more complex for modelling [2]. On the other hand, decision engineering applied to
chemical engineering leads us to define decision parameters (wj, pj, qj, Vj, Vj) in the form of
intervals or with fuzzy membership functions. In the case of a complex process, the natural
difficulty to define simple preferences leads to adapt decision engineering techniques. The
objective is to search for the decision capable of staying the best (or among the best) in spite of
small parameter variations. So, adaptation techniques and new projections in the decision aid
science will be needed to solve the problem of decision robustness.
It seems that our paper proves even the natural raison d'étre of the symbiosis between
scientific disciplines apparently very different, such as the science of decision making and the
chemical science in the area of controling various proc: . Now, chemical engineering has the
possibility to solve its multiobjective decision problems and decision engineering can adapt their
techniques to chemical industrial problems.
We hope by conviction that this present process control related research constitute a first
step towards the axiomatization of the multicriterion control of the real time decision-
making processes, whose industrial application is promising and innovative.
6 REFERENCES
[1] Massebeuf S., Fonteix C., Pla F., "Modélisation et optimisation multicritére du procédé de
polymérisation en émulsion du styréne", Pla-Schrauwen ed., /7@mes Journées Scientifiques du
Club Emulsion, 1998
[2] Derot B, J. Gareau, L.N. Kiss and J.-M. Martel: "The solver of Volvox Multicriterion Table",
in Multicriteria Analysis, J. Climaco(Ed.), pp. 113-126, Springer-Verlag, 1997
[3] C. FONTEIX, R. VIENNET, I. MARC :"A new experimental design for multicriteria
optimization: application to a biochemical synthesis". 2nd Int. Symposium on Mathematical
modelling and simulation in agriculture
and bio-industries, IMACS, 7-9 mai 1997, Budapest, (1997) 57-62
[4] Perrin E., Mandrille A., Oumoun M., Fonteix C., Marc I.: "Optimisation Globale par
Stratégie d’Evolution : Technique Utilisant la Génétique des Individus Diploides". RAIRO-
Recherche Opérationnelle, 31, 2, pp. 161-201, 1997
[5] Viennet R., Fonteix C., Mare L: "Multicriteria Optimization Using a Genetic Algorithm for
Determining a Set Pareto". Int. J. of Systems Science, 27, 2, pp.255-260, 1996
[6] N. KUEHM, C. FONTEIX, I. MARC : "Fuzzy dynamic experimental sests for the
determining of the validity of a state estimation".(Détermination de la validité d'une prédiction
d'état par Plans d'Expériences Dynamiques Flous.) RAJRO-APII-JESA Journal Européen des
Systémes Automatisés, 30, 5 (1996) 737-753
[7] Fonteix C., Bicking F., Perrin E., Marc I: "Haploid and Diploid Algorithms, a New
Approach for Global Optimization: Compared Performances". Int. J. of Systems Science, 26, 10,
pp. 1919-1933., 1995
[8] Viennet R., Fonteix C., Marc L: "New Multicriteria Optimization Method Based on the Use
of a Diploid Genetic Algorithm: Example of an Industrial Problem". Lecture Notes on Computer
Science. Artificial Evolution: European Conference, AE 95, Brest, France, Springer-Verlag Ed.
1063, pp. 120-127, 1995
[9] Courcoux P., Qannari E. M., Melcion J. P., Morat J.: "Optimisation Multiréponse :
Application 4 un Procédé de Granulation d’Aliments". Récents Progrés en Génie des Procédés :
Stratégie Expérimentale et Procédés Biotechnologiques, 36, 9, pp. 41-47, 1995.
[10] Horn J., Nafpliotis N., Goldberg D.E., "A niched Pareto genetic algorithm for
multiobjective optimization", Proc. of the first Conference on Evolutionary Computation, IEEE
World Congress on Computational Intelligence, vol. 1, pp. 82-87, 1994
[11] Rousseau A., Martel J.-M. : "Environnemental Assessment of on Electric Transmission
Line Project: A MCDM method", in Applying Multiple Criteria Aid for Decision to
Environmental Management, Paruccini (Ed.), Kluwer Academic Publishers, pp. 163-165, 1994
[12] Roy B., et Bouyssou D., "Aide multicritére 4 la décision: méthodes et cas", Economica,
Paris, 1993
[13] Fonseca C.M., Fleming P.J., "Genetic algorithms for multiobjective optimization:
formulation, discussion and generalization", Proc. of the fifth Int. Conf. On Genetic Algorithms,
Morgan-Kauffman, pp. 416-423, 1993
[14] Goldberg D.E., "Genetic algorithms in search, optimization, and machine learning",
Addison-Wesley, 1989
[15] Logothetis N., Haigh A., "Characterizing and optimizing multi-response processes by the
Taguchi method", Quality and Reliability Engineering International, 4, pp.159-169, 1988
[16] Schaffer J.D., "Multiple objective optimization with vector evaluated genetic algorithms",
In Grefenstette J.J. editor, Proc. Int. Conf. On Genetic Algorithms and their Applications, pp.
93-100, 1985
[17] Brans J.P., Marechal B. and Vincke Ph., "PROMETHEE: A new family of outranking
methods in MCDM", JFORS'84, North Holland, pp.477- 490, 1984
[18] Derringer G., Suich R., "Simultaneous optimization of several response variables", Journal
of Quality Technology, 12, pp. 214-219, 1980