Scholarly article on topic 'Minimizing Dimensional Variation and Robot Traveling Time in Welding Stations'

Minimizing Dimensional Variation and Robot Traveling Time in Welding Stations Academic research paper on "Materials engineering"

CC BY-NC-ND
0
0
Share paper
Academic journal
Procedia CIRP
OECD Field of science
Keywords
{Assembly / joning / "spot welding" / robots / "dimensional quality" / "cycle time" / optimization / "branch and bound"}

Abstract of research paper on Materials engineering, author of scientific article — Johan S. Carlson, Domenico Spensieri, Kristina Wärmefjord, Johan Segeborn, Rikard Söderberg

Abstract Complex assembled products as an automotive car body consist of about 300 sheet metal parts joined by up to 4000 spot welds. In the body factory, there are several hundred robots organized into lines of welding stations. The distribution of welds between robots and the welding sequences have a significant influence on both dimensional quality and throughput. Therefore, this paper proposes a novel method for quality and throughput optimization based on a systematic search algorithm which exploits properties of the welding process. It uses approximated lower bounds to speed up the search and to estimate the quality of the solution. The method is successfully tested on reference assemblies, including detailed fixtures, welding robots and guns.

Academic research paper on topic "Minimizing Dimensional Variation and Robot Traveling Time in Welding Stations"

Available online at www.sciencedirect.com

ScienceDirect

Procedia CIRP 23 (2014) 77 - 82

Conferenceon Assembly Technologies and Systems

Minimizing Dimensional Variation and Robot Traveling Time in Welding

Stations

Johan S. Carlsona, Domenico Spensieria*, Kristina Warmefjordb, Johan Segebornc, Rikard

Soderbergb

a Geometry and Motion Planning, Fraunhofer-Chalmers Centre, Chalmers Science Park, SE-412 88, Gothenburg, Sweden

b Geometry Assurance, Wingquist Laboratory at Chalmers, SE-412 96, Gothenburg, Sweden c Volvo Car Corporation, Gothenburg Sweden * Corresponding author. Domenico Spensieri, Tel.: +46 31 772 4252; fax: +46 31 772 4260. E-mail address: domenico.spensieri@fcc.chalmers.se

Abstract

Complex assembled products as an automotive car body consist of about 300 sheet metal parts joined by up to 4000 spot welds. In the body factory, there are several hundred robots organized into lines of welding stations. The distribution of welds between robots and the welding sequences have a significant influence on both dimensional quality and throughput. Therefore, this paper proposes a novel method for quality and throughput optimization based on a systematic search algorithm which exploits properties of the welding process. It uses approximated lower bounds to speed up the search and to estimate the quality of the solution. The method is successfully tested on reference assemblies, including detailed fixtures, welding robots and guns.

© 2014 Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).

Selectionand peer-reviewunderresponsibilityof thelnternational Scientific Committee of 5th CATS 2014 in the person of the Conference Chair Prof. Dr. Matthias Putz matthias.putz@iwu.fraunhofer.de

Keywords: Assembly, joning, spot welding, robots, dimensional quality, cycle time, optimization, branch and bound

1. Introduction

Complex assembled products as an automotive car body consist of about 300 sheet metal parts joined by up to 4000 spot welds. In the body factory, there is typically welding lines with several hundred robots organized into geometry welding stations and re-spot stations. The distribution of welds between robots and the welding sequences have a significant influence on both dimensional quality and throughput. Previous research by Segeborn et al. [1] has shown that math based automatic line balancing algorithms and methods for distributing the welds in re-spot stations, can improve cycle time up to 25%. In re-spot stations the geometry is considered to be independent of welding sequence. Furthermore, Genetic Algorithms (GAs) have proven to be an effective and efficient

approach for multi-criteria (simultaneous) optimization of welding sequences with respect to dimensional variation and robot traveling time, [2,3]. However, the GA approach cannot easily be extended to multi-robot stations, welding configuration alternatives, minimizing of robot coordination loses, estimation of the gap to the optimal sequence. These are all important aspect for optimization of geometry stations that are already part of the algorithms for re-spot stations. Therefore, this paper proposes a novel method based on a systematic search algorithm which exploits specific properties of the process behind, namely the welding operation. This informed algorithm allows creating approximate lower bounds in order to speed up the search and to estimate the quality of the solution obtained, in terms of gap to the optimum. Since the evaluations of these solutions are time consuming, the

2212-8271 © 2014 Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/3.0/).

Selection and peer-review under responsibility of the International Scientific Committee of 5th CATS 2014 in the person of the Conference

Chair Prof. Dr. Matthias Putz matthias.putz@iwu.fraunhofer.de

doi:10.1016/j.procir.2014.03.199

algorithm also strives to minimize the number of such computations.

Improving product quality and equipment utilization have both economic and environmental impact and thereby supports sustainable manufacturing. A more detailed discussion on sustainability measures and awareness in production used at Chalmers University of Technology can be found in Johansson et al. [4].

This paper has the following outline; In Section 2 the effect on geometrical variation of welding sequence is discussed and a novel algorithm for optimization is presented. Section 3 focuses on collision free path planning and optimization of robot traveling time. In Section 4, the algorithms are combined in a new way to optimize both quality and traveling time simultaneously. The optimization algorithm is tested on two reference assemblies in Section 5 and finally the conclusion are drawn in Section 6.

2. Geometrical Variation and Welding Simulation

In this Section we describe variation simulation of assembled products joined with discrete welding points. A novel method for optimizing welding sequence with respect to quality is presented.

2.1 Variation simulation model

Variation simulation aims at predicting variation and offsets in critical dimensions of assembled products or subassemblies. The inputs to such a simulation are, in the case of non-rigid simulation, part meshes, material properties, locating schemes, tolerances of the parts and fixtures, and information about the joining process - i.e. joining method, position of the joining points and joining sequence.

2.2 Locating schemes and assembly fixtures

Locating schemes are used to lock parts and assemblies theoretically in space. The theoretical locating points are in reality represented by locating features on the parts and the fixtures, i.e. holes, slots, pins and planes. For rigid analysis, six locating points are used to lock six degrees of freedom. For non-rigid assemblies, additional support points are used to represent clamps and supports in physical fixtures (Söderberg et al. [5]). Figure 1 shows a rigid location scheme, a non-rigid locating scheme and a welding fixture.

Fig. 1. Rigid and non-rigid location schemes and a welding fixture.

Since parts and fixture elements are afflicted with variation, resulting in positioning errors and gaps between the parts to be joined, a sub-assembly will not be nominal after joining.

Therefore, variation needs to be considered in joining simulation. Also the sequence in which they are joined can have a significant influence on the final geometrical quality, see Figure 2 (Warmefjord et al. [6]). The figure shows geometrical variation in a number of inspection points for an A-pillar assembly.

Fig. 2. The geometrical effect of different welding sequences in a number of inspection points.

2.3 Variation simulation including contact modelling and joining sequence

The variation simulation used in this paper is based on Monte Carlo simulation enclosed in the software RD&T, RD&T Technology, [7]. Here, a total sensitivity matrix is implicitly defined by a FEA-based simulation model describing all mating conditions, kinematic relations and nonrigid behaviour. Direct Monte Carlo simulation combined with FEA is a standard technique for variation prediction of compliant parts. However, since a large number of replications are required for achieving satisfactory accuracy, the method is very time-consuming if a new FEA calculation is to be executed in each replication. Liu and Hu [8] presented a technique called Method of Influence Coefficients (MIC) to overcome this drawback, and this method is used in this paper.

To avoid having the parts virtually penetrate each other during assembly, a point-based contact modelling algorithm is used (Warmefjord et al. [6]). The joining sequence is also taken into consideration. For a non-rigid assembly, the following steps are necessary to take into consideration when predicting variation and deviation in critical dimensions of the final assembly:

■ Step 1: The parts are positioned/clamped in their fixtures, and over-constrained locating systems (i.e. clamps) are applied. Forces are applied to clamp nonnominal parts.

■ Step 2: The parts are joined together in a pre-defined joining sequence. The gaps in the joining points are closed, one by one.

■ Step 3: After the last joint is set, the assembly is unclamped and is allowed to spring back.

More comprehensive explanations of these steps can be found in Warmefjord et al. [9, 10].

2.4 Optimization of welding sequence

Evaluating welding sequences is even with the state of the art simulation method described above very time consuming. Therefore, in this section we present a novel method for welding sequence optimization which tries to keep the evaluations to a minimum. For moderate cases as with 7 welds the combinatorial possibilities are 5040. The basic idea of the method is the experimental evidence that simultaneous welding gives better quality outcome than welding one point after another, see Figure 2 and Warmefjord et al. [9]. This evidence is supported by our experimental data in many cases. Such a fact can be used as an assumption to provide "almost exact" lower bounds for the geometrical quality of a welding sequence. The possibility to create these lower bounds naturally suggests the use of branch and bound-based (B&B) algorithms to optimize variation in welding sequences. In Figure 3, a B&B tree is depicted, illustrating a small example with 4 welding points, where the sequences are drawn in black and the lower bound value in grey: the sequence 3 1, for example, means that first point 3 is welded, than point 1, thereafter points 2 and 4 are welded simultaneously. One value representing geometrical variation may be proportional to standard deviation, e.g. 6a, or other measures can be used such as Mean Squared Error, MSE.

Fig. 3. The geometrical effect of different welding sequences in a number of inspection points.

We adopt, here, an eager approach: lower bounds are computed only when B&B nodes are expanded, in order to decrease the number of simulations needed. Nodes expansion is usually done by depth-first strategies since they perform well for B&B, anyway we use here a best-first queuing system because we do not want to explore the entire tree. Note that best-first is a suitable approach regarding memory requirements, since in these applications the number of points is low. In this small case, 9 sequences out of 24 were evaluated to find the optimal solution. Generally for a case with N welds the algorithm can in best case find the optimal solution in O(N2) evaluations with an eager approach, and O(N) with a lazy one. We will come back to the efficiency and effectiveness of the proposed algorithm on the industrial reference cases in Section 5.

3. Path Planning and Robot Traveling Time

To program robot motions and find a time efficient sequence for a spot weld robot station like the one shown in Figure 4 is a time consuming and really challenging task. The programming can be done by using the real physical robot or

by using a computer model of the robot. The former is called teach in programming, where an operator uses a joystick to control the motions of the robot in order to specify the robot program. The latter is a digitalized version of teach in programming and is called off-line programming (OLP). OLP has many advantages such as better equipment utilization and faster introduction of new product models. However, one of the main advantages is that by using OLP one can further reduce the programming time and improve the program by math-based algorithms for motion planning and combinatorial optimization.

The overall approach for optimization of travel time in welding robot station consists of two major steps, (i) task planning to find promising configurations that can reach each weld and (ii) sequence optimization and motion planning to select one solution for each weld and connect them together by efficient motions and in a sequence minimizing traveling time. However, before this is described in more detail (Section 3.4) we introduce automatic collision free path planning and exact solving of generalized traveling sales problems (GTSP).

3.1 Automatic Path Planning

Automatic path planning addresses the problem of finding collision free motions of moving objects in cluttered environments. Complete algorithms are of little industrial relevance due to the complexity of the problem (PSPACE-hard for polyhedral models (Canny, [11]). Instead, sampling based techniques trading completeness for speed and simplicity are the methods of choice. Common for these methods are the needs for efficient collision detection, nearest neighbor searching, graph searching and graph representation. The two most popular methods are; Probabilistic Roadmap Methods (PRM) and Rapidly-Exploring Random Trees (RRT). The PRM samples randomly among all configurations of an object, keeps the collision-free samples, and then connects those pair wise if the straight line path between them is collision-free (Bohlin and Kavraki, [12]). The RRT incrementally builds two trees from the start and the goal configurations respectively. In each step an attractor is generated at random and trees are expanded from the nearest node towards the attractor. The iteration stops when the trees overlap (LaValle and Kuffner, [13]). Inspired by both these probabilistic methods Fraunhofer-Chalmers Centre (FCC) has developed a deterministic path planner that adaptively adjusts a grid in the configuration space. However, the methods proposed in this paper can be based on any efficient path planner algorithm.

3.2 GTSP solved by Dynamic Programming

The traveling salesman problem (TSP) can be solved exactly by easy-to-implement algorithms, which however have high computational complexity. One of these is due to Held and Karp, see [14] and is based on dynamic programming: it improves the complexity of brute force search, from O(n!), to O(n2 2n). The Bellman equations describing the solution are:

Y 0) = Ci;

KSJ) = minieSX{jî{7(S\{/},i) + ctJ}

(Eq. 1)

Here Cij is the cost to move between nodes i and j, whereas the function y{S,j) represents the minimum cost to move from node 1, visiting all nodes in S and ending at node j. It is immediate to note that its practical use is limited to small size problems: on a modern PC it would run on instances of about 20 nodes. The same approach can be used to solve exactly the GTSP. Only a straightforward generalization is necessary, leading to the following recursive formulas:

Y(j)[m] = min!eiv(0 c1;[l,m] Y(S,j)[m] = minieSX{ji[min!eW(i){KS\{7},i)[Z] +cy[l,m]}}

Here Cy[l, m] represents the cost to move from the /h node in group i to the mth node in group j. At the same strength as in Eq. 1, y(5,7)[m] refers to the minimum cost of visiting group 1, all groups in S and ending at the mth node in group j, and N(i) is the set of all nodes contained in group i. For small problems (of the same order as for the TSP, i.e. circa 20), the above equations allow to solve to optimality even some variants of the GTSP, like the ones fixing home and/or end group, asymmetric instances, and fixing the order of some parts of the tour. The latter can be used within the B&B framework to combine bounds for geometrical variation and traveling time for different welding sequences, as will be described in Section 4.

3.3 Task planning

Due to the spin symmetry normal to the surface of the welding points there exists many ways to weld each point without affecting the quality. Also, different inverse robot solutions exist as elbow up or down. These alternative solutions play an important role in traveling time optimization. Therefore, in the task planning step we calculate all collision free alternatives with 5 degrees resolution on the spin angle to reach each weld, see Figure 4.

Fig 4. Task planning: Highlighting 2 configurations out of 32 possible.

3.4 Path Planning and Sequence optimization

The next step is sequencing where it is decided in which order and with what alternative configuration the robot should perform the spot weld to minimize the cycle time. This is a generalization of the classical Traveling Salesman Problems (GTSP) since it is a grouped problem with node alternatives. A straight-forward algorithm would be to calculate, by path

planning, the full cost tensor {mtjk[] between alternative j and I of spot i and k and then apply a GTSP solver to find the optimal solution. However, the path planning operations between pairs of welds are much more time consuming than solving the GTSP itself. Therefore, we propose a lazy method minimizing the number of calls to the path planner. First a lower bound cost matrix rni]ki = ||9tj — 8kl || „ is calculated based on linear motions between the weld configurations 8tj and 8ki in the six dimensional configuration space for the robot. By running a GTSP algorithm on the cost matrix m a route R will be generated. In the next step, the cost matrix is updated by path planning of all pair wise configuration included in the route R. By iteratively running these steps, the optimal solution can be found with minimal calls to the path planner. The iterative sequence optimization is illustrated on a schematic case in Figure 5.

Fig. 5: A four group generalized TSP problem with obstacles. Left figure shows the optimal routing before path planning. Right figure shows the rerouting after a number of iteration considering collisions with the obstacles (Mark et al. [15]).

4. Minimizing Variation and Traveling Time

One way to simultaneously measure geometrical quality, q, and cycle time, t, of a given sequence of N welding points is by their weighted sum. By changing the weight, different solutions on the Pareto front can be found. Limitations and extension of weighting methods for multi criterias can be found in e.g. Athan and Papalambros [16]. Since q and t are different physical entities, we normalize them by using two reference values: the best lower bound for geometrical quality, Q0, and the best lower bound for cycle time, T0. Note that these values are obtained at the root node of the B&B tree, where no part of the welding sequence is assigned. The objective function, therefore becomes

f = af + 0.-a)±

where a £ [0,1] gives the possibility to differently weight the quality measure q for that specific welding sequence and its correspond cycle time t. Note that q > Q0 and t> T0. This objective function can be used along the B&B. At each node i of the tree, representing a subsequence of M points (P;i,..., P;M), we have

fi = «+ (1 ->l Qo 'To1

where qt is the quality measure obtained by sequentially welding (Pjl, ...,P;M), and then simultaneously welding the remaining N-M points, and t; is the cycle time optimizing a GTSP where the first M points have the fixed order {Pi1, ■■■, Pim). The former is computed by variation simulation software, whereas the latter can be optimally computed by small modifications to the algorithm described in Section 3.4. In Figure 6 a graphical representation of a fictive B&B node consisting of 5 welding points is illustrated. Two points P1 and P3 are fixed, whereas P2, P4 and P5 can be taken in any order, after P1 and P3. The arrows represent the sequential constraints between points just described. Furthermore, each point can be welded by several robot configurations alternative, depicted by the black dots within each circle in the Figure. In Figure 7 a possible solution is illustrated.

HOME ,-.

Fig. 6. The sequential constraints for the B&B node P1Pi.

Fig. 7. A solution to the problem modeled in fig. 6. 5. Industrial Reference Assemblies and Results

In this section the proposed methods for welding sequence optimization are evaluated on two reference assemblies. Reference Assembly I is an A-pillar and Reference Assembly II is a part of the front floor, see Figures 8 and 9. A more comprehensive description can be found in Segeborn et al. [3]. Both assemblies consist of two sheet metal parts which are joined by seven weld points. For Reference Assembly I, a number of 30 instances of Part I and Part II have been physically measured on single part level, using 52 inspection points on part I and 17 inspection points on Part II. Then, using the RD&T variation simulation tool, all 30x30=900 assemblies combining Part I and Part II are virtually built. For Reference Case II, the variation in each contact point and locator is statistically simulated and a number of 3000 assemblies are built. The output variation q is the dimensional variation of the resulting assembly, defined as the root mean square of the 6a values of all nodes of the assembly mesh. For both assemblies all 5040 sequences are simulated and there quality are plotted in Figure 10.

Fig. 8. Industrial Reference Assembly I. The assembly is virtually measured using the fixture of Part I, including only support S3.

Reference Assembly N

Fig. 9. Industrial Reference Assembly II. The assembly is virtually measured using the fixture of Part I.

Fig. 10. Geometrical variation (6sigma) for all 5040 welding sequences of Reference Assembly I and II.

The algorithm described in Section 2.4 has been applied to the two assemblies to calculate a welding sequence maximizing the global geometrical quality. Different settings for the B&B stop criterion have been tested and the results are reported in the Table 1.

#evaluations Quality % above optimum Rank (out of 5040)

100 1,739 2,2 48

200 1,710 0,4 4

46 1,739 2,2 48

Table 1. Geometrical quality optimization by B&B applied to Reference Assembly I.

The first two rows refer to runs of the algorithm stopping after 100, respectively 200 iterations and saving the best sequence obtained. In the third row, the resulting sequence corresponds to the first expanded leaf in the B&B tree, reached after 46 expanded nodes. Note that this sequence is the same as the one stopping after 100 iterations. The quality improvement obtained by running the algorithm 200 iterations, against 100,

is clear. However, it is an open discussion whether it is worth the computational effort, and it will most probably depend from case to case. The second reference case, whose results are reported in Table 2, shows a larger improvement in quality, from running 200 against 100 iterations, with respect to the previous test case.

Acknowledgements

This work was carried out within the Wingquist Laboratory VINN Excellence Centre, supported by the Swedish Governmental Agency for Innovation Systems (VINNOVA). This work is also part of the Sustainable Production Initiative and the Production Area of Advance at Chalmers University of Technology.

#evaluations Quality % above optimum Rank (out of 5040)

100 1,398 5,5 405

200 1,348 1,7 100

29 1,411 6,5 523

Table 2. Geometrical quality optimization by B&B applied to Reference Assembly II.

Since the method for optimizing quality works well, using few evaluations finding top quality sequences close to optimum solution on the two test cases, the next step is to test the multi-criteria algorithm described in Section 4. In Figure 11 an overview of the algorithms building blocks and their connections are illustrated. The optimization core consists of the 3 blocks: B&B, Variation Simulation (VS), and AGTSP Solver (GTSP). When a solution is obtained from that, path planning is performed, exploiting the lazy philosophy as in Section 3.4. Data are updated and a new iteration is run with the new information gathered. Note that when the VS block is not present, i.e. a=0, the algorithm becomes as the one in Section 3.4, whereas, cutting the GTSP block leads to optimizing with a=1. In the latter case, the groups sequence is fixed at any iteration: only the choice of the nodes within the groups is left to optimization.

Fig. 11. General framework with building blocks.

The results for the simultaneous optimization of geometrical quality and cycle time for three different values of a (the max number of evaluations is set to 100) are presented in Table 3. Simulations are stopped after 3h computing and the tradeoff between quality and cycle time is evident.

a Geometrical quality Cycle time (s)

1.0 1.739 7.672

0.9 1.811 4.944

0.0 1.859 4.441

Table 3. Tradeoff between quality and cycle time for Assembly 1.

6. Conclusion and Future Research

An important first step to optimize geometry stations with respect to both dimensional quality and cycle time has been proposed and successfully applied to reference assemblies. Next step will be to extend to full geometry stations making it necessary to include load balancing and coordination of several robot sharing space and the welding work.

References

[1] Segeborn, J., Segerdahl, D., Ekstedt, F., Carlson, J. S., Mikael Andersson, Söderberg, R., "An Industrially Validated Method for Weld Load Balancing in Multi Station Sheet Metal Assembly Lines", ASME Journal of Manufacturing Science and Engineering, 2013.

[2] Xie, L. Shelley, and Ching Hsieh. "Clamping and welding sequence optimisation for minimising cycle time and assembly deformation.", International Journal of Materials and Product Technology 17.5: 389399, 2002.

[3] Segeborn, J., Carlson, J.S., Wärmefjord, K., Söderberg, R.,"Evaluating Genetic Algorithms on Welding Sequence Optimization with Respect to Dimensional Variation and Cycle Time", ASME 2011 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, August 28 - September 31, Washington, DC, USA, 2011.

[4] B. Johansson, A. Dagman, E. Rex, T. Nyström, M. Knutson Wedel, J. Stahre, and R. Söderberg, "Sustainable production research: awareness, measures and development," International Journal of Sustainable Development, vol. 4, no. 11, pp. 95-104, Sep. 2012.

[5] Söderberg, R., Lindkvist, L. and Dahlström, S., "Computer Aided Robustness Analysis for Compliant (non-rigid) Assemblies", Journal of Engineering Design, Vol. 17, No. 5, pp. 411-428, 2006.

[6] Wärmefjord, K.,, Söderberg, R. and Lindkvist, L.. Tolerance simulation of compliant sheet metal assemblies using automatic node-based contact detection. in Proceedings of IMECE2008, Boston, USA.

[7] RD&T Technology "RD&T Reference Manual", Mölndal, Sweden

[8] Liu, S.C. and Hu, S.J., "Variation Simulation for Deformable Sheet Metal Assemblies Using Finite Element Methods", Journal of Manufacturing Science and Engineering, Vol. 119, No. 3, pp. 368-374, 2007.

[9] Wärmefjord, K.,, Söderberg, R. and Lindkvist, L., "Strategies for Optimization of Spot Welding Sequence With Respect to Geometrical Variation in Sheet Metal Assemblies", in Proceedings of ASME International Mechanical Engineering Congress & Exposition, Vancouver, Canada, 2010.

[10] Wärmefjord, K., Söderberg, R. and Lindkvist, L., Variation Simulation of Spot Welding Sequence for Sheet Metal Assemblies. in Proceedings of NordDesign2010 International Conference on Methods and Tools for Product and Production Development, August 25-27, Gothenburg, Sweden. 2010.

[11] Canny, J.F., The Complexity of Robot Motion Planning, MIT Press, 1988.

[12] Bohlin R. and Kavraki, L.E., Path Planning Using Lazy PRM, In: Proc. IEEE Int. Conf. On Robotics and Automation, 2000.

[13] LaValle, S.M., and Kuffner, J.J., Randomized Kinodynamic Planning. In: Proc. IEEE Int. Conf.on Robotics and Automation, 1999.

[14] Held, M., and Karp, R.M., "A dynamic programming approach to sequencing problems", in Proceedings of the 1961 16th ACM National Meeting, pages 71.201-71.204, New York, NY, USA.

[15] A. Mark, A., Bohlin, R., Segerdahl, D., Edelvik, F., Carlson, J.S., "Optimization of Robotized Sealing Stations in Paint Shops by Process Simulation and Automatic Path Planning"International Journal of Manufacturing Research, 2013.

[16] Athan, T. W., & Papalambros, P. Y., "A note on weighted criteria methods for compromise solutions in multi-objective optimization", Engineering Optimization, 27(2), 155-176, 1996.