Discrete representation of non-dominated sets in multi-objective linear programming

Discrete representation of non-dominated sets in multi-objective linear programming

 

Lizhen Shao a , ∗, Matthias Ehrgott b

  • Key Laboratory of Advanced Control of Iron and Steel Process (Ministry of Education), School of Automation and Electrical Engineering, University of Science and Technology Beijing, Beijing 10 0 083, China
  • Department of Management Science, Lancaster University Management School, Bailrigg, Lancaster LA1 4YX, United Kingdom
a r t i c l e i n f o

Article history:

Received 3 September 2015

Accepted 2 May 2016

Available online 27 May 2016

Keywords:

Multi-objective optimisation

Linear programming

Non-dominated set

Discrete representation

a b s t r a c t

In this paper we address the problem of representing the continuous but non-convex set of non- dominated points of a multi-objective linear programme by a finite subset of such points. We prove that a related decision problem is NP-complete. Moreover, we illustrate the drawbacks of the known global shooting, normal boundary intersection and normal constraint methods concerning the coverage error and uniformity level of the representation by examples. We propose a method which combines the global shooting and normal boundary intersection methods. By doing so, we overcome their limi- tations, but preserve their advantages. We prove that our method computes a set of evenly distributed non-dominated points for which the coverage error and the uniformity level can be guaranteed. We ap- ply this method to an optimisation problem in radiation therapy and present illustrative results for some clinical cases. Finally, we present numerical results on randomly generated examples.

© 2016 Elsevier B.V. All rights reserved.

 

1.Introduction

Multi-objective linear programming (MOLP) problems arise in many real world applications of Operations Research. Due to the presence of multiple conflicting objective functions, there is in gen- eral no single optimal value but an infinite continuous set of non- dominated objective function vectors. A non-dominated objective function vector (or non-dominated point) is the image Cx of an ef- ficient solution x of the MOLP. An efficient solution x is a feasible solution of the MOLP such that there is no other feasible solution which is at least as good as xinallobjectivesandstrictlybetterinatleastone.

In practical applications, a decision maker has to choose a most preferred point from the set of all non-dominated points. Since the non-dominated set of an MOLP with p objectives is in gen- eral a set of dimension p − 1 , this is a difficult task. This is the motivation for computing a finite representative subset of non- dominated points: selecting from a finite set is easier and there is a wide variety of decision analysis methods for this task, see, e.g., Greco, Ehrgott, and Figueira (2016).Naturally,therepresen-

More information on the data used in this paper can be found at doi: 10.17635/ lancaster/researchdata/76 .

∗ Corresponding author. Tel.: +86 1062332926.

E-mail addresses: lshao@ustb.edu.cn (L. Shao), m.ehrgott@lancaster.ac.uk (M. Ehrgott).

http://dx.doi.org/10.1016/j.ejor.2016.05.001 0377-2217/© 2016 Elsevier B.V. All rights reserved.

tative set should be chosen in such a way that its cardinality is not too big, that each non-dominated point is close to at least one representative point (coverage error) and that the representative points are uniformly distributed and not too close to each other (uniformity level). We will give formal definitions (originally due to Sayin, 20 0 0 ) in Section 2.

In this paper, which is an extension of an earlier conference pa- per ( Shao & Ehrgott, 2007 ), we present a method to compute a finite representative subset of the non-dominated set of an MOLP. In Section 2 we formally define multi-objective linear programmes, review the definitions of the criteria for a good representation and point out a relationship to so called dispersion problems of loca- tion theory. We formulate the finite representation problem as a multi-objective optimisation problem over the non-dominated set of an MOLP and use the relationship to dispersion problems to prove that a decision version of our formulation is NP-complete. Section 3 is devoted to a review of the literature on methods to compute finite representations, with focus on the global shoot- ing method ( Benson & Sayin, 1997 ) , the normal boundary inter- section method ( Das & Dennis, 1998 ) and the normal constraint method ( Messac, Ismail-Yahaya, & Mattson, 2003 ). We provide ex- amples that show that they all may fail the coverage property or provide bounds on the uniformity level of the representation. Our own method, combining features of the global shooting and normal boundary intersection methods is described in Section 4 . In Section 5 we prove bounds on the coverage error and the uniformity level of the representative set computed by our method under certain assumptions on the underlying MOLP. We also show that the method works in polynomial time in the bi-objective case. Section 6 provides some examples, presents an application of our method to an MOLP arising in radiotherapy treatment planning and provides numerical results on randomly generated examples. The numerical results provide empirical evidence for the theoreti- cal results of the previous section. Section 7summarisesourcon-tributionsandidentifiessomequestionsforfurtherresearch.

2.Multi-objectivelinearprogrammingandfiniterepresentationofaset

We shall use the following notation for the comparison of vec- tors y 1 and y 2 ∈ R p . We write y 1  y 2 if yfor all k = 1 , . . . , p.

We use y 1y 2 to indicate y 1  y 2 but  , whereas y 1 < y 2 means y.Weshallalsousethenotation

 

Amultipleobjectivelinearprogrammingproblemincompactmatrixnotationcanbewrittenastheoptimisationproblem

min {Cx  : xX } ,                                                                                      (1)

where C ∈ R p× n is a p × n matrix of objective function coeffi- cients, the rows c k T , k = 1 , 2 , . . . , p of which are the coefficients of p linear functions c k T x , k = 1 , 2 , . . . , p. We assume that X ⊂ R n is a nonempty compact polyhedral set of feasible solutions in decision space R n . The feasible set Y in objective space R p isdefinedby

Y = {Cx  : xX } .                                                                                        (2)

Since the image of a convex polyhedron under a linear map is also a convex polyhedron (see, e.g., Rockafellar, 1970 ), it follows that Y is also a nonempty compact convex polyhedron. In this paper, we shall make the further assumption that Y N contains at least two distinct points. Otherwise the non-dominated set is a singleton, and the MOLP (1)istrivial.

Wenextdefineefficientsolutionsandnon-dominatedpoints.

Definition 2.1. Feasible solution x 0X is a (weakly) efficient so- lution to MOLP (1) , if there does not exist any xX such that Cx ( < ) ≤ Cx 0 . The set of all efficient solutions of MOLP (1) will be denoted by X ( W ) E and is called the (weakly) efficient set in de- cision space. Correspondingly, y 0 = C x 0 is called a (weakly) non- dominated point and Y (W ) N = {Cx  : xX E } is the (weakly) non- dominated set in objective space for problem (1).

Theorem 2.2isfundamentalinmulti-objectivelinearprogram-ming.

Theorem 2.2 ( Isermann, 1974 ) . Feasible solution x 0X is an effi- cient solution of MOLP (1) if and only if there exists a λ ∈ R >p  such that

T 0                         T Cx                                                                                         (3)

for all xX.

We recall that it is well known (see, e.g., Naccache, 1978; Yu & Zeleny, 1975 ) that the non-dominated set Y N of an MOLP (1) is connected and consists of the union of faces of Y . It will be con- venient to refer to maximal faces of Y that belong to Y Nasnon-dominatedfaces.

Definition 2.3. Let FY be a face of Y such that FY N and F is maximal with this property with respect to set inclusion. Then Fiscallednon-dominatedface.

Weshallnowintroducethenotionofrepresentationofaset,thequalityofsucharepresentation,andthecomplexityofcom-puting one.

Let Z ⊂ R p be a set and let RZ be a finite subset. R is then called a discrete representation of Z . We are interested in “how well” R represents Z . Sayin (20 0 0) defines coverage, uniformity, and cardinality as the three attributes of quality of discrete represen- tations. According to these three attributes, a good representation needs to contain a reasonable number of points, should not miss large portions of Z,andshouldnotcontainpointsthatareveryclosetoeachother.Thisleadstothefollowingdefinitions.

Definition 2.4. Let   0 be a real number and d be a metric. Rep- resentation R of Z is called an -representation of Z if for every zZ , there exists some rR such that d ( z , r )  .

Definition 2.5. Let δ  0 be a real number and d be a metric. Rep- resentation R of Z is called a δ-uniform representation of Zif

.

Wenotehere,thatthroughoutthepaperwewillalwaysusetheEuclideandistanceasmetric.

Sayin (20 0 0) proposes measures to quantify the three quality attributes. The number of points contained in a representation ob- viously measures its cardinality. The coverage error  and unifor- mity level δ aredefinedasfollows.

The coverage error  signifies how precisely the set Z is being represented by the discrete representative subset R . It can be writ- ten as:  = max min d(z, r) . zZ rR

Thus, how well a fixed zZ is covered is determined by the closest point to z in the representation R . For the entire set Z , the coverage error depends on how well an arbitrary element of Z is covered. Therefore, the coverage error  is equal to the maximum coverage error for individual points in Z.

The coverage error is closely related to what we informally term the coverage propertyofamethodtogenerateadiscreterepresen-tation,namelytheproperty,thateverynon-dominatedpointcanpotentiallybeselectedasarepresentativepoint.

The uniformity level of representation R can be measured by the distance between a pair of closest distinct points of R.Thusitcanbeexpressedas

1 2 .

 

Foradiscreterepresentation,asmallnumberofpoints,lowcoverageerror,andhighuniformitylevelaredesirable.

From the discussion above, we can formally write the problem of finding a discrete representation of Zasamulti-objectiveopti-misationproblem

min | R |                                                                                                      (4)

max                 min d(r i , r j )                                                                      (5)

 

min max min d(z, r)    (6) zZ rR

s.t. RZ, | R | < ∞ .                                                                                (7)

Multi-objective optimisation problem (4) –(7) is the problem of locating a small number of points in set Z such that the points are far apart from each other, yet each point of Z is close to at least one point of R . Such problems are known as dispersion problems in location theory. In the context of this paper, therefore, Z would be Y N,thenon-dominatedsetofanMOLP.

This relationship between dispersion problems and multi- objective decision making was in fact first observed by White (1991) : “… in the context of multiple objective problems, it may be that it is not only the maximal dispersion which is of interest, but also the representativeness of the solution obtained for the set Zasawhole.” However,heconsiderstheproblem

max { min{  d(x, y ) : x, yR } : RZ, | R | = p} .

Thus, White (1991) considers only the uniformity measure, (5) , keeps the cardinality constant and ignores the coverage error. He also restricts Z to a finite set. More relevant is the work by Baur and Fekete (2001).Theyconsider,amongstothers,thesocalledpuredispersionproblem

max { min{  d(v , w ) : v , wR } : RZ, | R | = p} , where Z is a polygonal region in R 2 and d(v , w ) is the geodetic distance between v and w . They prove Theorem 2.6.

Theorem 2.6. Unless P = NP, there is no polynomial-time approxi- mation scheme for pure geometric dispersion.

From this result it follows that the decision problem: Given a polygonal region Z , an integer p , and a constant δ, does there exist a subset R of Z with cardinality p such that min{  d(v , w ) : v , wR }  δ? is NP-complete. Hence the canonical decision problem of the discrete representation problem stated in (4) –(7) , given Z , an in- teger p and nonnegative scalars  and δ, does there exist a finite subset R of Y N of cardinality p such that the coverage error is at most  and the uniformity level is at least δ is NP-complete, since it is NP-complete for the specific choice of  such that d ( z 1 , z 2 )   for all z 1 , z 2Z.

To relate this to the problem of finding discrete representa- tions of the non-dominated set of an MOLP, we notice that, as mentioned in Section 1 the non-dominated set of an MOLP as de- fined in (1) is the union of finitely many non-dominated faces of Y . Therefore, finding discrete representations of Y Nisequivalenttofindingdiscreterepresentationsofaunionofpolytopes.

To complement this hardness result, we notice that for discrete bi-objective optimisation problems Vaz, Paquete, Fonseca, Klam- roth, and Stiglmayr (2015) provide algorithms that solve the dis- crete representation problem for known Y N in time polynomial in | Y N | and | R | for various combinations of coverage error, unifor- mity level, and -indicator as quality measures. Eusébio, Figueira, and Ehrgott (2014) provide algorithms to compute δ-uniform or -representations for bi-objective integer network flow problems, but they do not analyse their complexity. In Section 4 we prove that our method computes a discrete representation of the non- dominated set of a bi-objective linear programme by solving O (| R |) linear programmes without the knowledge of Y N.

3.Asurveyofexistingmethods

3.1. Brief literature review

We can categorise the methods for finding discrete represen- tations of the non-dominated set of multi-objective optimisation problems into two groups, one is based on the knowledge of X E and the other works without the knowledge of X E.

Based on the knowledge of X E , Sayın (2003) proposes a proce- dure to find discrete representations with specified coverage errors. The procedure also specifies the uniformity level of the representa- tions. Knowledge of X E can, however, not be assumed when solving a multi-objective optimisation problem. Therefore most methods work without the knowledge of X E . A recent survey of methods for computing discrete representations in multi-objective optimisation can be found in Faulkenberg and Wiecek (2010).Inthefollowing,wefocusonmethodsthatarerelatedtotheconceptsoftheglobalshootingandthenormalboundaryintersectionmethods.

Benson and Sayin (1997) propose a global shooting method to find a representation of the non-dominated set. This method has the coverage property, but it cannot guarantee the uniformity, i.e. a bound on the uniformity level, of the representations it gener- ates. We will discuss the global shooting method in more detail in Section 3.2 . Das and Dennis (1998) propose a normal boundary in- tersection (NBI) method for finding non-dominated points for gen- eral (nonlinear) multi-objective optimisation problems. It uses the convex hull of the individual minima (CHIM) of the p objectives as a reference plane. Evenly distributed reference points are placed on the CHIM and for each reference point a non-dominated point is computed by solving a single objective optimisation subproblem. While the method generates evenly distributed non-dominated points, some parts of the non-dominated set may be missed, i.e. it does not have the coverage property, a problem caused by the use of the CHIM. We will illustrate the method and its limitations in Section 3.3.

Based on the NBI and the -constraint methods, Ismail-Yahaya and Messac (2002) propose the normal constraint (NC) method. Instead of an equality constraint used in the subproblems of the NBI method, the NC method uses an inequality constraint to re- duce the feasible set in the objective space of the multi-objective optimisation problem. Later, Messac et al. (2003) propose the nor- malised normal constraint (NNC) method. The NNC method works in a normalised objective space. However, both the NC and NNC methods have the same drawback as the NBI method, because they use the CHIM as reference plane. Realising this limitation of using the CHIM, Messac and Mattson (2004)improvetheNNCmethodbyusinganextendedCHIMinsteadoftheCHIMasreferenceplane.However,uniformitylevelandcoverageerrorofthediscreterepre-sentationarenotmeasured.

Martinez, Sanchis, and Blasco (2007) improve the NNC method for bi-objective optimisation problems. The improvements are based on two types of techniques, i.e., nonlinear optimisation and genetic algorithms. Sanchis, Martínez, Blasco, and Salcedo (2008) propose a new alternative method, the enhanced nor- malised normal constraint (ENNC) method for multi-objective optimisation problems. They present the formulation of a new reference plane that improves the original normalised normal constraint method of Messac et al. (2003) using two approaches: a redefinition of the anchor points, i.e., vertices of the reference plane and an exact linear transformation between the objective space and the normalised objective space. Mueller-Gritschneder, Graeb, and Schlichtmann (2009) propose a successive approach to systematically build up the representative non-dominated set. The approach is based on the generation of so-called trade-off limits. Motta, Afonso, and Lyra (2012) propose a novel modified procedure which is similar to the approach of Mueller-Gritschneder et al. (2009) and they claim that their method when applied to the NBI and NC methods for more than two objectives overcomes some of their deficiencies. Logist and Van Impe (2012) give theoretical insights in the conditions under which the NBI and the ENNC are able to generate the same candidate non-dominated points. Hancock and Mattson (2013) propose the smart normal constraint method for generating a “smart” non-dominated set. However, with the idea of smart they do not consider uniformity any more. Another recent paper dedicated to generating equidistant repre- sentations in bi-objective optimisation proposed by Faulkenberg and Wiecek (2012)isbasedonthe-constraintmethod.

It is worth mentioning that there are many approximation methods for multi-objective optimisation (the reader is referred to Ruzika & Wiecek, 2005)thatcomputesomenon-dominatedpoints.However,sincetheirgoalistoapproximatethewholenon-dominatedsetandtheydonotaimatfindingevenlydistributednon-dominatedpointstheymayyieldbadrepresentationsintermsofcoverageerroranduniformitylevel.Suchmethodsarethereforenotmentionedintheabovesurvey.

We now present more details of some of the methods mentioned above. We note that although they are originally formulated

forgeneralmulti-objectiveoptimisationproblems,weuseMOLPnotationthroughoutforconsistency.

3.2. The global shooting method

As indicated in Section 1 , we assume that the feasible set Y in objective space is bounded. The global shooting method defines polyhedron Y: Cx  y  yˆ for some xX} , where yˆ is chosen as a point so that for all yY N it holds that y  yˆ.  For in- stance, yˆ can be chosen as the anti-ideal point y AI the components of which are defined by y AIk  := max{  y k , yY } , k = 1 , . . . , p. Benson (1998) has shown that Y   has dimension p and that  and Yhavethesamenon-dominatedset.

The global shooting method then constructs a simplex S that contains  . A subsimplex Sˆ  of S (the non-dominated set of S ) is chosen as the reference plane. Then, a discrete sample of reference points are placed on Sˆ  and the method “shoots” from yˆ towards each reference point as far as possible while remaining in   This is achieved by solving the single objective linear programme (the global shooting subproblem (8) ) . max t

s.t.                     yˆ + t(qyˆ )  Cx

(8) t                      0 x                   X,

where q is a reference point. Thus, a set of points on the bound- ary of Y   is computed and each reference point q i corresponds to a boundary point y i = yˆ + t i (qyˆ ) of Y   , where t i is the optimal value of linear programme (8) for q = q i . However, not every inter- section point y i is non-dominated. Therefore it needs to be checked whether the intersection point y i isdominatedornotbysolvingthefollowinglinearprogramme.

min λT y

s.t. y  y i (9) y  ,

where λ ∈ R p is an arbitrary strictly positive vector. It is most con- venient to choose λ = e, the vector of all ones. Point y i is non- dominated if and only if the optimal value of (9) is λT y i . In case y i is dominated, an optimal solution to (9) is non-dominated and is added to the representative set in place of y i .

Fig. 1 illustrates the global shooting method for an MOLP prob- lem where Yistheconvexhullofvertices(3,7),(2,9),(3,12),(6,11),(8,9),(9,7),(10,4)and(6,5).Thenon-dominatedextremepointsoftheMOLPare(2,9),(3,7),(6,5)and(10,4).Thereare13ref-erence pointsonthe referenceplane which isthelinesegment

Fig. 2. Unevenly distributed representative non-dominated points obtained by the global shooting method.

from ( −2,12) to (10,0). For each reference point, we solve prob- lem (8) and obtain 13 intersection points with  . Among the 13 intersection points, nine are non-dominated and four are dom- inated. The latter are shown as crosses in Fig. 1 . After solving (9) , the nine non-dominated intersection points are kept, whereas solving (9) for the four weakly non-dominated intersection points yields non-dominated point y = (2 , 9).Thustherepresentativesetsconsistoftheninenon-dominatedintersectionpointsandnon-dominatedpoint(2,9).

For the MOLP case, which is the topic of this paper, the global shooting method simply involves solving two linear programmes for each reference point. It has the coverage property because it puts reference points on Sˆ  and Y NSˆ  + R p  . We refer to Benson and Sayin (1997)forproofsofthepropertiesoftheglobalshootingmethod.

Next,weprovideaseriesofMOLPsthatshowsthattheuni-formitylevelofthediscreterepresentativesetcomputedbytheglobalshootingmethodcanbearbitrarilysmall,eveniftherefer-encepointsareequidistantlydistributedonthereferenceplane.

Example 3.1. Consider the MOLP problem min{ Cx  : Ax  b, x ∈ R n }with

 

 

C  ,

 

where M  1.

Fig. 2 shows the feasible set Y in objective space for M =

9 . The non-dominated set is the line segment from point (M − 1 , M + 1) to point ( M , 1). We use the global shooting method to obtain a set of non-dominated points. M + 2 reference points which are evenly distributed on the reference plane (the line segment between (0 , M + 1) and (M + 1 , 0))areused.Theyare

(0 , M + 1) , (1, M ), (2 , M − 1) , . . . , (M − 1 , 2) , ( M , 1) and (M + 1 , 0).

We obtain M +1correspondingnon-dominatedpointsshownin

Fig. 2 . They are . ,

and ( M,1).TheEuclideandistancebe-

tween− the “first” two−      non-dominated points (M − 1 , M + 1)and

is d 1 = √ 22  1andthedistancebe-

M +                                   2 2 M1

tween the “last” two non-dominated        points ( 3 M3  M−1 , M 23 + M4M1 1 )

(M

andtimes ( M d , 1  .1) As is  Md 2  approaches= +31 M)M1   2infinity, +1 so that d 1  dapproaches 2 is equal zero, while d 2

approaches infinity. Therefore, as M → ∞ , the uniformity level tends to 0 and the coverage error tends to infinity. This clearly

showsthattheglobalshootingmethodcannotguaranteeanypos-itiveuniformityleveloranyboundedcoverageerror.

3.3. The Normal Boundary Intersection (NBI) method

We first assume the individual minima of the functions c k T x over X are attained at x k for k = 1 , 2 , . . . , p. Let y k = Cx k and let y I  , c Tp  x p )T  denote the ideal point. We note that in the original paper, Das and Dennis (1998) consider a multi- objective optimisation problem in which Y is translated by y I sothattheidealpointbecomes0.Tobeconsistentwiththedescrip-tionoftheothermethodsinthispaper,wedonotworkwiththetranslatedproblem.

The convex hull of the individual minima (CHIM) is then de- fined as the set of all convex combinations of the individual global minima of the objective functions, i.e. conv {Cx  k : k = 1 , . . . , p} . A set of evenly distributed reference points on the CHIM is generated and for each of them a NBI subproblem is solved to find the far- thest point on the boundary of Y along the normal nˆ of the CHIM pointing towards the ideal point. The NBI subproblem for a given reference point q is the linear programme max t

s.t.                    q + t nˆ = Cx

(10)

t  xX.

Fig. 3 shows how the NBI method works for the same MOLP example with two objectives as in Fig. 1 . For this example, all the points obtained are non-dominated and the non-dominated set is uniformly covered. For problems with more than two objectives, Das and Dennis (1998) state that the method may overlook por- tions of the non-dominated set. They claim that these overlooked areas are likely near the periphery of the non-dominated set. How- ever, in Example 3.2 we provide an MOLP with p =3objectivesthatshowsthatthenormaltotheCHIMmaynotalwaysbeposi-tiveanddoesthereforenotfindanynon-dominatedpointonthelargestpart(anyfacet)ofthenon-dominatedset.Infact,inthisexample,theNBImethodonlyfindsnon-dominatedpointsthatbe-longtotheCHIM.Inotherwords,theNBImethoddoesnothavethecoverageproperty.

Although Das and Dennis (1998)claimthattheNBImethoddoescomputeevenlydistributednon-dominatedpoints,theydonotprovideboundsontheuniformitylevelofthecomputedrep-resentativeset.

3.4. The Normal Constraint (NC) and related methods

Like the NBI method, the normal constraint method ( Ismail-

Yahaya & Messac, 2002)alsousestheCHIM(calledtheutopiaplaneinthatpaper) asthe reference plane toputreferencepoints.

Fig. 5. The normalised objective space of the NNC method.

However,foreachreferencepoint,NCsolvesadifferentsubprob-lemtocomputeacorrespondingnon-dominatedpoint.Thesub-problemisdescribedasfollows:

min y p

s.t. N k T (y  , for all k = 1 , . . . , p − 1 (11) y   Y,

where N k = y py k and y k = Cx k is a vector such that y kk  is the minimal value of the kthobjective.

The inequality constraint N k T (yq )  0reducesthefeasiblesetintheobjectivespaceoftheMOLP.Therefore,thesubproblemisactuallyminimisingthelastsingleobjectiveinthereducedfeasiblesetinobjectivespacesubjecttothereducedobjectivespace.

Fig. 4 shows how the NC method works for the same MOLP example with two objectives as in Fig. 3.

For handling disparate objective magnitudes, scales, or ranges of objective function values for the different objective functions, Messac et al. (2003) propose the normalised normal constraint (NNC) method. The NNC method is the same as the NC method, except that the subproblems are formulated in a normalised ob- jective space, see Fig. 5.

Similar to the NBI method, the NC and NNC methods do not provide a guarantee that the generated set of non-dominated points will represent the non-dominated set well if p  3. Realising this limitation, Messac and Mattson (2004) improve the NC method by using an extended CHIM instead of the CHIM as

Fig. 6. and the non-dominated set in Example 3.2 .

reference plane. They claim that this improvement will gener- ate evenly spread non-dominated points over the entire non- dominated set. However, as Example 3.2 shows, this is not necessarily always the case. Messac and Mattson (2004) also generalise the subproblem. Instead of minimising the last single objective, the subproblem can minimise any single objective. Messac and Mattson (2004) expect that the same result is ob- tained no matter which objective is minimised. However, this is not always true. Example 3.2canbeusedtoshowthatsolvingdifferentsubproblemgivesdifferentresults.

Example 3.2. We consider the linear relaxation of a binary assign- ment problem with three objectives. The cost coefficients of the three linear objective functions are shown as matrices c 1 , c 2 and c 3,

 

 

c ,

 

 

c

We define Y: Cx  y  yˆ for some xX} and choose yˆ as (21, 21, 21), which is greater than the anti-ideal point (20, 20, 20) in this example. Fig. 6 shows Y   astheconvexhullofthe16extremepoints(21,21,21),(11,21,21),(21,9,21),(21,21,10),(11,11,14),(11,11,21),(15,9,21),(11,21,14),(21,14,10),(19,21,10),(13,21,11),(19,14,10),(21,9,17),(21,11,14),(15,9,17),(13,16,11).

The four points y 1 = (11 , 11 , 14) , y 2 = (15 , 9 , 17) , y 3 =

(19 , 14 , 10) and (13, 16, 11) represented by circles are the non- dominated extreme points of  . The non-dominated set consists of a line segment from point y 1 to point y 2 and a triangular facet which is the convex hull of y 1 , y 3 and(13,16,11).

The three dots at y 1 , y 2 and y 3 in Fig. 6 are the unique non- dominated points at which the individual minima of the three ob- jectives are attained. Hence the CHIM is the convex hull of these three points. The normal of the CHIM is nˆ = (1 , −40 , −28) or nˆ = (−1 , 40 , 28) , which are not positive. Placing evenly distributed ref- erence points on the CHIM, the results of the NBI, NC and NNC methods are shown in Fig. 7 , 8 and 9 , respectively. In those fig- ures, circles represent the reference points, while dots represent non-dominated points generated by the respective methods. It is worth noting here, that neither Das and Dennis (1998) nor Messac et al. (2003)specify proceduresforchecking thenon-dominance

Fig. 7. Reference points and non-dominated points generated by the NBI method in Example 3.2 .

Fig. 8. Reference points and non-dominated points generated by the NC method in Example 3.2 .

Fig. 9. Reference points and non-dominated points generated by the NNC method in Example 3.2 .

of points generated by solving subproblems (10) and (11) , respec- tively. However, in this example, the optimal solutions of both sub- problems for reference points located on the CHIM are those ref- erence points themselves. For testing the methods, we did, there- fore, implement the same non-dominance check that we use for the RNBI method described in Section 4below.

We did not obtain any non-dominated points in the interior of the non-dominated facet by any of the three methods. For the NBI method, this occurs because the angle between the non-dominated facet and the CHIM is approximately 151.5 > 90 degrees. But even for the NC and NNC method, the subproblems only gener- ate non-dominated points on the (relative) boundary of the non- dominated facet and on the non-dominated edge between y 1 and y 2 .Hence,allthreemethodsclearlymissthemajorportionofthenon-dominatedset.

WealsoappliedtheimprovedNCandNNCmethodswithex-tendedCHIM.Wechosetheextendedreferenceplanebigenoughsothattheprojectionofthenon-dominatedsetiscontainedintheplane.FortheNCmethod,wealsosolveddifferentsubproblemsfor the same setofreferencepointsto seeiftheyyieldthesameresultor not,theydidnot.We omitdetailedresultsandfigures,butnotethatallgeneratednon-dominatedpointsareontheboundaryofthenon-dominatedset.

While Benson and Sayin (1997) have shown that the global shooting method has the coverage property, Example 3.2showsthattheNBIandNCmethodsdonothavethisproperty.Thisre-sultindicates,thatingeneralthesemethodswillnotbeabletocompute-representationsforarbitrary.

4.Therevisednormalboundaryintersectionmethod

As discussed in Section 3.2 above, the global shooting method satisfies the coverage property, and the NBI method ( Section 3.3)cangenerateevenlydistributednon-dominatedpoints.Therefore,weproposearevisedNBImethodwhichcombinesthesetwoap-proaches.Therevisednormalboundaryintersection(RNBI)methodprovidesaprioriguaranteesonboththeuniformitylevelandcov-erageerror.

On the one hand, instead of the CHIM, the RNBI method uses the non-dominated subsimplex Sˆ  of the simplex S that is used in the global shooting method ( Benson & Sayin, 1997 ) as the refer- ence plane on which to put equidistant reference points. By doing this, it retains the coverage property. On the other hand, by solving single objective subproblems similar to (10),theRNBImethodgen-eratesasetofnon-dominatedpointswithaknownpositiveboundontheuniformitylevel.

The RNBI method involves choosing a reference plane, placing equidistant reference points on the plane and computing the in- tersection points of rays emanating from the reference points in a direction normal to the reference plane with the feasible set Yinobjectivespace.Finally,thenon-dominanceofthecomputedin-tersectionpointsischeckedbysolvingalinearprogramme.Inthefollowingparagraphsweexplaineachofthesestepsindetail.

Constructing the reference plane: The reference plane is the same as in the global shooting method of Benson and Sayin (1997).Forcompleteness,wedescribethedetailsofthisconstructionhere.First,let

β := min{  e T y : yY } ,                                                                        (12)

where e ∈ R p isavectorinwhicheachentryis1.

Next, we define p + 1 points v k ∈ R p , for k = 0 , 1 , . . . , p. Let v 0 = y AI and, for k = 1 , 2 , . . . , p and l = 1 , 2 , . . . , p,let

v lk   y AIl  ,                      , ifif  llkk,,                                             (13)

As Benson and Sayin (1997) have shown, the convex hull Sof

{ v k : k = 0 , 1 , . . . , p} is a p -dimensional simplex, and S contains Y.

The subsimplex of S given by the convex hull Sˆ  of { v k : k = 1 , 2 , . . . , p} is the non-dominated set of S . It is contained in the hyperplane { y ∈ R P : e T y = β} , where β∗ is the optimal value of (12) . This hyperplane supports Y and Y Natalloptimalsolutionsof

(12).

Placing equidistant points on the reference plane: Clearly, for p = 2 , Sˆ  is a line segment. Proposition 4.1 shows that in the general case of p > 2 objectives, Sˆ  is a p − 1dimensionalsimplexwithequaledgelength.

Proposition 4.1. The reference plane Sˆ  is an equilateral simplex.

Proof. Let us consider two vertices v i and v j , i, j ∈ { 1 , . . . , p} and assume without loss of generality that i < j.

v iv j = (v i1  , . . . , v ii  , . . . , v ij  , . . . , v ip  ) T − (v 1j  , . . . , v ij  , . . . , v jj  , . . . , v pj  ) T

− −

(y AI1  , . . . ,y AIi  , . . . ,β + y AIj  − e T v 0 ,. . . , y AIp  ) T

Fig. 10. Equidistant reference points on the reference plane.

Fig. 11. The revised normal boundary intersection method.

 y AIi  , . . . , y AIj

−(β + y AIj  e T v 0 ) , 0 , . . . , 0) T

= (0 , . . . , β e T v 0 , . . . , −β + e T v 0 , . . . , 0) .

Thus the distance√             between any two vertices of Sˆ  isequalto

d(v i , v j ) = 2 (β − e T v 0 ) .

For p = 3 , Fig. 10 shows how a triangular lattice can be used to generate equidistant points on Sˆ .  In the general case of p objec- tives, the i th reference point q i isgivenby

p  

q i  ki  v k k =1

ifrom0to1with

where k  and k = 1 finitek  set . By of varying equidistant αk  points on the a fixed increment of ηk a reference plane can be generated. For the three objective case in Fig. 10 ηk = 0 . 25 . We denote dsthedistancebetweentwoclosestreferencepoints.

Computing the intersection points: Given a reference point q on Sˆ ,  the RNBI subproblem searches for the point in Y that is clos- est to the reference point along the normal direction e . This is achieved by solving the RNBI subproblem (14) . min t

s.t.           q + teY     (14) t

Fig. 11 illustrates the RNBI method for the same MOLP example that was used in Figs. 1 and 3 to illustrate the global shooting and the NBI method. The reference plane and the placement of reference points are the same as for the global shooting method. However the search for intersection points with Y uses a direction normal to the reference plane, as in the NBI method. Since due to our assumption of boundedness of Y (and therefore S ), the linear programme (14) cannot be unbounded, it is either infeasible or has an optimal solution. In fact, there are three scenarios for the solu- tion of (14) , see also Fig. 11.

  1. (14) is This occurs if and only if there is no intersec- tion between the ray { q + te : t 0 } and Y.
  2. (14) has an optimal solution t ∗ if and only if the ray { q + te : t 0 } and Y In this case, y = q + te belongs to the boundary of Y.Twosubcasesmayoccur.
  • The intersection point y ∗ is
  • The intersection point y ∗ isnon-dominated.

Because in the case that (14) has an optimal solution, the iden- tified intersection point y ∗ may or may not be non-dominated it is necessary to check the status of y ∗. A simple non-dominance fil- ter can be applied to the set of all generated intersection points to eliminate some of them, see Messac et al. (2003) . This method has the advantage of being fast, but it may accept some dominated points (often ones close to Y N)asnon-dominated.

An exact method to test the non-dominance of an intersection point y ∗ is similar to (9) and is provided by Proposition 4.2 , a proof of which can be found in Ehrgott (2005).

Proposition 4.2. Let y Y and let   Then y Y N if and only if yis an optimal solution to the linear programme (15) ,

min λT y

s.t.                (15) yY.

5.PropertiesoftheRNBImethod

InthissectionweinvestigatethepropertiesoftheRNBImethod,inparticularweproveboundsontheuniformitylevelandthecoverageerroroftherepresentativesetgeneratedbytheRNBImethod.

First of all, we notice that for every non-dominated point yˆ of the MOLP (1) there is a point qˆ ∈ Sˆ  such that (14) has an optimal solution tˆ  such that qˆ + tˆ e  = yˆ.  This point is simply the projection in direction e of yˆ onto Sˆ .  Hence it follows, similar to Benson and Sayin (1997)fortheglobalshootingmethod,thattheRNBImethodhasthecoverageproperty.Nextweshallinvestigatetheuniformitylevel.

Let F be a non-dominated facet of Y with normal nˆ.  Since the normal of Sˆ  is e by definition, the angle between Sˆ  and Fisgivenby

=       e T nˆ                            nˆ 1 + ··· + nˆ p

cos θ        =                                                                .                                      (16)

Since, and Definition 2.3thenormal

nˆ to F is an element of R >p  we have

1 + ··· +

.                                                                                                                     (17)

Therefore

(18)

and θ isintherangeof

For p 4 and for p , i.e., the range of angles between the reference plane Sˆ  and non-dominated facets increases with p.

Fig. 12. Analysis of distances for a 2D example.

Now for the equidistant reference points with distance ds that the RNBI method places on Sˆ  it follows that the corresponding non-dominated points in the representative set have a distance of √

d s/ cos θ < pd  s.

Fig. 12 shows an example with two objectives ( p = 2 ) so that non-dominated facets are line segments. F 1 is a non-dominated facet, while F 2 is a weakly non-dominated facet. Angles between non-dominated facets and the reference plane must be smaller than  , whereas the angle between the reference plane and the weakly non-dominated facet F 2 is . The distance between two closest non-dominated√  points obtained by the RNBI method is between ds and 2ds  .

Since for an MOLP (1) with p objectives we have shown that for the distance d between two closest non-dominated points in the representative set it holds that ds  d < √ pds , we have   ds . Thus, Theorem 5.1isproved.

Theorem 5.1. Let R be the representative subset of Y N obtained by the RNBI method. Then R is a ds-uniform representation of Y N .

The width w (B ) of a convex set B ⊂ R p is defined as the small- est Euclidean distance between two supporting hyperplanes of B . Since this means that the width of any convex set of dimension less than p is 0 but we want to measure the width of projections of subsets of Y N onto Sˆ ,  we also define the width of a convex set on a hyperplane Sˆ ,  w Sˆ  (B ) as the minimal distance between two parallel supporting hyperplanes perpendicular to Sˆ .

It is well known ( Yu & Zeleny, 1975 ) that the non-dominated set of MOLP (1)isthefiniteunionofmaximalnon-dominated

facesto the of reference Y , i.e. Y  Nplane = ∩ Kk  = Sˆ1 .   FY k N .p   Now,can  thenlet  Ybe Np   bewritten the  projectionas the union of  Yof N

K polytopes  O k on Sˆ ,  where O k is the projection of F k onto Sˆ .  Using this notion, we can now prove the bound of √ pds  on the coverage error of the discrete representation generated by the RNBI method under an assumption on the width of sets O k.

Theorem 5.2. Let ds be the distance between reference points and assume that each set O k in the projection Y Np  of Y N satisfies w Sˆ  (O kds, then the representative set R obtained by the revised NBI method

is a        pds -representation of Y N .

Proof. Let Q be the set of reference points for which solving (14) results in the set of representative points R . Since the width of O k , k = 1 , . . . , K on the reference plane Sˆ  is greater than or equal to ds , it follows that QO k for all k = 1 , . . . , K.

Let yY N and let o be the projection of y to Sˆ .  Clearly oO k for some k ∈ { 1 , . . . , K} . Therefore, there must exist qQO k such that d ds because the distance of reference points is ds . Clearly, there also exists rR which is the representative point generated by q . From the analysis of distances between represen- tative points, it follows that d.

Therefore, for any non-dominated point yY N,theremustexist

rR such that d.

We remark that the proof of Theorem 5.2 also shows that, even in the case that not all O k have width at least ds , any non- dominated facet of Y N that satisfies the condition w Sˆ  (O kds will be present in the representation. Hence it implies that only non- dominated faces that are “too small” because either their dimen- sion is smaller than p − 1 or their width is smaller than √pds  can  

bemissed.

Theorems 5.1 and 5.2 quantify the quality of the discrete rep- resentation generated by the RNBI method in terms of coverage error and uniformity level. They also clarify that the only param- eter of the method is the distance ds between two closest refer- ence points. Since Sˆ  is an equilateral simplex, this then determines the number of reference points and hence an (as we shall see very weak) upper bound on the cardinality of the discrete representa- tion, the lower bound ds on the uniformity level, and the upper √

bound pds  on the coverage error. By increasing ds , the cardinality will increase, the uniformity level will decrease and the coverage error will decrease. We emphasise that the RNBI method does not allow independent control of the three measures of quality of the discrete representation. Nevertheless, it is to the best of our knowl- edge the first method that allows the computation of a discrete representation of the non-dominated sets of MOLPs with guaran- teed coverage error and uniformity level, and without the knowl- edge of Y N.

Finally, we consider the special case of MOLPs with p = 2 ob- jectives. Let y lex (1, 2) and y lex (2, 1) denote the two lexicographically optimal non-dominated points and o lex (1, 2) and o lex (2, 1) their pro- jections on the reference plane. Then for any reference point q such that q 1 < olex 1  (1 , 2) or q 1 > olex 1  (2 , 1) we have that (14) is either infea- sible or its optimal solution yields a dominated point. On the other hand, all reference points with olex 1  (1 , 2)  qolex 1  (2 , 1) yield non- dominated points. Hence, to compute a representation R of cardi- nality | R | = r one places r equally spaced reference points on the line segment from o lex (1, 2) to o lex (2, 1) and solve (14) for these refer- ence√  points. The resulting representation will have coverage error 2δ and uniformity level δ,where

δ = d(olex  (1 , 2) , olex  (2 , 1) ) / (r − 1) .

Theorem 5.3. For any MOLP with p = 2 objectives,√  such that Y is bounded, the RNBI method computes a δ-uniform-representation of cardinality r for Y N in time polynomial in the size of the MOLP and linear in r.

6.Examplesandnumericalresults

In this section, we first apply the RNBI method to Examples 3.1 and 3.2,whichweusedbeforetoillustratethattheglobalshootingandNBImethodscannotprovideaguaranteeontheuniformitylevelorthecoverageerror.Then,wepresentanapplicationoftheRNBImethodtothebeamintensityoptimisationproblemarisingintheplanningofradiotherapytreatmentofcancer.Finally,we present theresultsofnumericaltests ofthe

Fig. 13. The non-dominated points generated by the RNBI method for Example 3.1 .

Fig. 14. The set of non-dominated points generated by the RNBI method for Example 3.2 .

RNBImethodonrandomlygeneratedlinearprogrammeswithbetweenthreeandeightobjectives.

In Example 3.1 we apply the RNBI method with the same set of reference points as we used for the global shooting method. The RNBI method generates six non-dominated points as shown in Fig. 13 . These six points are evenly distributed on the non- dominated set and the distance between the two closest points in the representation is 1.9079.√  The guaranteed uniformity level in this example is of course 2.  According to Theorems 5.1 and 5.2 , the√  distance between√ √ two  closest representative points is at least √ 2 and at most 2 · 2 = 2 .HencetheRNBImethodcomputesa2-uniform2-representationofthenon-dominatedset.

Next, we apply the RNBI √method  to Example 3.2 . We use ref- erence points with distance 2 as in the previous example. This results in 325 reference points. The RNBI method generates 33 in- tersection points, 10 of which are non-dominated points. These are shown in Fig. 14 . We note that all non-dominated points generated by the method are on the non-dominated facet defined by vertices (11, 11, 14), (19, 14, 10) and (13, 16, 11). No non-dominated point on the edge between (11, 11, 14) and (15, 9, 17) is generated. This is of course expected due to Theorem 5.2,becausethewidthoftheprojectionofthenon-dominatededgeontothereferenceplaneisequalto0.Thus,theRNBImethodcorrectlycomputesarepresen-tationofthenon-dominatedfacet,withoutguaranteeinganycov-erageofthenon-dominatededge.

We shall now apply the RNBI method to a multi-objective linear programme derived from the so-called beam intensity optimisation problem in the planning of radiation therapy treatment for cancer. This serves to illustrate the potential of the RNBI method in prac- tical application. We refer the reader to Ehrgott, Güler, Hamacher, and Shao (2009) for more details on optimisation problems in this domain. The goals of radiation therapy are to deliver a uniform dose as close as possible to a prescribed dose to the planning target volume (which contains the tumour to be treated), while at the same time sparing the surrounding normal tissues and organs at risk as much as possible from the harmful effects of radiation. Given the number of beams and beam directions from which to irradiate the patient, beam intensity optimisation computes beam

8 12 40 60 α 20 20 10 0 −10 −20 −30 −40 β α β α β

Fig. 15. The non-dominated points generated for the AC, PR and PL cases (from left to right) when 153 reference points are used.

Fig. 16. The non-dominated points generated for the AC, PR and PL cases (from left to right) when 378 reference points are used.

Table 1

The number of reference points, intersection points with Y , non-dominated points, distance be- tween closest reference points and CPU time in seconds for the application of the RNBI method to an MOLP in radiotherapy treatment planning.

  RP IP NP ds CPU
AC 378 92 92 0 .89 27 .562
PR 378 148 117 4 .64 31 .437
PL 378 154 137 3 .17 110 .736
AC 153 42 42 1 .36 12 .351
PR 153 62 48 7 .09 20 .684
PL 153 62 56 4 .84 52 .153

intensity profiles that yield the best dose distribution for treating a particular patient under consideration of clinical and physical constraints. Due to the conflicting goals, the beam intensity opti- misation problem can, for example, be formulated as an MOLP in the form provided in Shao and Ehrgott (2008) . In this MOLP the objectives are to minimise the maximum deviations α, β, γ ofde-livereddosefromprescribedlowerboundsonthedosedeliveredtothetumourandfromprescribedupperboundsonthedosede-liveredtotheorgansatriskandothernormaltissue,respectively.

We use three clinical cases, an acoustic neuroma (AC), a prostate tumour (PR) and a pancreatic lesion (PL). These are or- dered according to the number of constraints and the number of variables. The RNBI method was implemented in Matlab 7.3 using the Matlab optimisation toolbox as LP solver and all tests were run on a PC with 2.5 gigahetrz processor speed and 4.0 gigabyte RAM. In Table 1 , we list the number of reference points (RP), the num- ber of intersection points generated (IP), as well as the number of non-dominated points (NP). The fourth column shows the distance dsbetweentwoclosestreferencepoints.Finally,theCPUtimeinsecondsisshown.Foreachcase,weusedtwodifferentnumbersofreferencepoints.

We notice that for all three cases, more than half of the refer- ence points did not generate intersection points. However, check- ing infeasibility of (14) was not time consuming in this application, but this may be an issue for other MOLP problems. Moreover, Table 1 shows that for the PR and PL cases not every intersection point is non-dominated. Hence, checking non-dominance by solv- ing linear programme (15)isnecessary,butmaybemoretimeconsumingthancheckinginfeasibility.Clearly,themaineffortinapplyingtheRNBI methodisinsolvingLPs.Hence,the timetakendependsonthenumberofreferencepoints,sinceatleastone(butatmosttwo)LPsneedtobesolvedforeachreferencepoint.

We show the non-dominated points of the three clinical cases in Figs. 15 and 16 . The generated non-dominated points clearly ap- pear to be evenly distributed. We also note that the distance ds of reference points and the distance of non-dominated points is mea- sured in Gray, the unit of radiation dose, because all three objec- tives are measured in this unit. This unit has a clear meaning for radiotherapy planners. Since each non-dominated point represents a potential treatment plan, their distance can be interpreted as a measure of the difference of the plans. Radiotherapy planners may have a good idea what dose value would constitute recognisable significant difference in plans. Hence they should be able to set a desired value for the distance d of non-dominated points,√ which  then implies that ds should be set to a value between d/ 3 and d to √guarantee  non-dominated points with a distance between d and 3d . The value of dtheninturndeterminesthenumberofreferencepointstobeused.

To conclude this section, we present numerical results for ran- domly generated MOLP problems. The examples are generated as follows. First, l points x i ∈ R p , i = 1 , . . . , l with x ran- domly distributed between 0 and 1 and x ip

  − 1) 2 are generated. Next, the convex hull of the lpoints

1 x i ,. , l is constructed. Now, we let the convex hull be the feasible set X in decision space, i.e., we use all facet defining in- equalities of the convex hull as constraints of the MOLP. Finally, we set the objective matrix C of the MOLP to be the identity ma- trix I . This procedure is employed in order to generate MOLPs with p -dimensional feasible set in objective space for which the non- dominated set has non-dominated facets. Since the randomly gen- erated points are distributed on the lower-left part of the sphere x n = (x 1 − 1) 2 + . . . + (x n −1 − 1) 2 , this does happen most of the time (always in our experiment). If, instead, we were to gener- ate the coefficients of the constraint matrix, right hand side vector and objective function matrix of (1) randomly, it would quite often be the case that Y and Y N have lower dimensions, in which case Theorem 5.2tellsusthattheRNBImethodmaynotwork.

Ten examples for selected combinations of l (number of points) and p (number of objectives, which is equal to the number n of variables) were solved for two different numbers of reference points (RP). The average number of constraints maswellastheresults,consistingoftheaveragenumbersofintersectionpoints(IP), non-dominatedpoints(NP),distancebetweentwoclosest ref-

Table 2

The number of objectives p , variables n , random points l , constraints m , reference points RP, intersection points IP on Y , non-dominated points NP, the distance ds between two closest reference points and the CPU time in seconds.

p = n l m RP IP NP ds CPU
3 30 56 78 42 .7 40 .4 0 .1851 0 .7811
      153 74 .8 70 .9 0 .1364 1 .3024
4 40 195 286 82 .4 74 .2 0 .2613 2 .1546
      816 202 .8 180 .9 0 .1926 5 .9825
5 50 748 727 85 .3 71 .4 0 .3511 7 .3551
      3080 293 .8 237 .5 0 .2587 25 .1582
6 60 3041 1378 83 .5 63 .0 0 .4339 16 .9257
      8723 355 .9 251 .8 0 .3197 131 .4911
7 70 12678 2003 57 .3 46 .0 0 .4988 93 .5675
      19292 350 .3 240 .8 0 .3675 810 .3551
8 80 53239 2281 41 .1 28 .7 0 .5687 445 .2241
      34122 259 .4 146 .7 0 .4191 6140 .5817

Fig. 17. The generated non-dominated points and feasible sets in objective space for two of the randomly generated examples with p = 3 .

erence points ds and the CPU time in seconds are listed in Table 2 . We also show the generated non-dominated points and the feasi- ble sets in objective space for two of the randomly generated ex- amples with (p, l, m ) = (3 , 30 , 56) in Fig. 17 . The one on the left has (RP,IP,NP, ds ) =(153,81,73,0.1273)andtheoneontheright(78,40,37,0.2122).

Table 2 shows that the CPU time clearly increases with the size of the MOLP (number of variables and constraints) and also the number of reference points, but as is usual in multi-objective op- timisation, the most significant impact on CPU time is the num- ber of objective functions p . Nevertheless, the RNBI method can be applied to MOLPs with up to eight objective functions. Under the assumption of Theorem 5.2 it is guaranteed to generate an evenly distributed representative set of non-dominated points. It is also apparent that as p increases, the number of reference points that do not generate an intersection point with Yincreasesdramatically.

7.Conclusionandfuturework

Inthispaper,wehaveaddressedtheproblemofgeneratingadiscreterepresentationofthenon-dominatedsetofamulti-objectivelinearprogramme.Wehaveproposedtherevisednor-malboundaryintersectionmethod,which,bycombiningfeaturesofthenormalboundaryintersectionmethodandtheglobalshoot-ingmethodovercomesdrawbacksofbothofthesemethodsandgeneratesanevenlydistributedsetofrepresentativepoints.Infact,wehavebeenabletoproveaprioriboundsontheuniformitylevelandthecoverageerrorofthediscreterepresentationundertheassumptionthatthenon-dominatedfacesofthefeasiblesetinobjectivespaceare“bigenough” relativetothedistanceofrefer-encepoints.TheRNBImethodisthefirstmethodforwhichqualityguaranteesintermsofbothcoverageerroranduniformitylevelhavebeenproved.Moreover,numericalresultsonintensityopti-misationproblems fromradiotherapytreatmentplanning andonrandomlygeneratedexamplesempiricallyconfirmthatthepointsofthediscreterepresentationareindeedevenlydistributed,andthatitisapplicableforMOLPswithuptoeightobjectives.

We have indicated in Section 6 above, that one of the problem- atic issues with the RNBI method is that many reference points lead to infeasible RNBI subproblems. Hence methods to decide whether (14) is feasible before solving it, would help eliminate un- necessary computation. This is a topic we will further investigate in the future. Moreover, we have seen that (14) may yield dom- inated intersection points with Y . In the case of p = 2 objectives this can easily be avoided by replacing the anti-ideal point y AI by the nadir point y N in the construction of the reference plane. This makes it easy to determine which RNBI subproblems will yield a non-dominated point and the RNBI method becomes an algorithm that runs in time linear in the cardinality of the representation. However, the determination of the nadir point is itself a difficult problem for p  3, so other techniques to replace the anti-ideal point by a better upper bound on Y Nareworthinvestigating.

In order to address more general problems, we will first gener- alise the method to MOLPs without the compactness assumption on Y , i.e., where Y may be unbounded, with Y Neitherboundedornot.Anotherinterestingextensiontobeaddressedinfuturere-searchisavariantoftheRNBImethodforconvexmulti-objectiveoptimisationproblems.

Acknowledgements

This work has been partially supported by Beijing Natu- ral Science Foundation [Grant number 4152034 ]; National Natu- ral Science Foundation of China [Grant numbers 810 0 0650 and 61370132]andNationalHigh-techResearchDevelopmentProgramofChina(863Program)[Grantnumber2013AA040705].

References

Baur, C. , & Fekete, S. (2001). Approximation of geometric dispersion problems. Algo- rithmica, 30 (3), 451–470 .

Benson, H. P. (1998). An outer approximation algorithm for generating all efficient extreme points in the outcome set of a multiple objective linear programming problem. Journal of Global Optimization, 13 (1), 1–24 .

Benson, H. P. , & Sayin, S. (1997). Towards finding global representations of the ef- ficient set in multiple objective mathematical programming. Naval Research Lo- gistics, 44 (1), 47–67 .

Das, I. , & Dennis, J. E. (1998). Normal-boundary intersection: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems.

SIAM Journal on Optimization, 8 (3), 631–657 .

Ehrgott, M. (2005). Multicriteria optimization (2nd). Berlin, Heidelberg: Springer .

Ehrgott, M., Güler, C., Hamacher, H. W., & Shao, L. (2009). Mathematical optimiza- tion in intensity modulated radiation therapy. Annals of Operations Research, 175 (1), 309–365. doi: 10.1007/s10479- 009- 0659- 4 .

Eusébio, A. , Figueira, J. , & Ehrgott, M. (2014). On finding representative non-domi- nated points for bi-objective integer network flow problems. Computers & Oper- ations Research, 48 , 1–10 .

Faulkenberg, S. L. , & Wiecek, M. M. (2010). On the quality of discrete representa- tions in multiple objective programming. Optimization and Engineering, 11 (3), 423–440 .

Faulkenberg, S. L. , & Wiecek, M. M. (2012). Generating equidistant representations in biobjective programming. Computational Optimization and Applications, 51 (3), 1173–1210 .

Greco, S., Ehrgott, M., & Figueira, J. (Eds.) (2016). Multiple criteria decision analysis: State of the art surveys. International series in operations research & management science (vol. 78) (2nd). New York: Springer. doi: 10.1007/b100605 .

Hancock, B. J. , & Mattson, C. A. (2013). The smart normal constraint method for directly generating a smart Pareto set. Structural and Multidisciplinary Optimiza- tion, 48 (4), 763–775 .

Isermann, H. (1974). Proper efficiency and the linear vector maximum problem. Op- erations Research, 22 , 189–191 .

Ismail-Yahaya, A. , & Messac, A. (2002). Effective generation of the Pareto frontier us- ing the normal constraint method. In Proceedings of the 40th aerospace sciences meeting and exhibit . Paper no. AIAA 2002-0178.

Logist, F. , & Van Impe, J. (2012). Novel insights for multi-objective optimisation in engineering using normal boundary instersection and (enhanced) normalised normal constraint. Structural and Multidisciplinary Optimization, 45 (3), 417–431 .

Martinez, M. , Sanchis, J. , & Blasco, X. (2007). Global and well-distributed Pareto frontier by modified normalized normal constraint methods for bicriterion problems. Structural and Multidisciplinary Optimization, 34 (3), 197–209 .

Messac, A. , Ismail-Yahaya, A. , & Mattson, C. A. (2003). The normalized constraint method for generating the Pareto frontier. Structural and Multidisciplinary Opti- mization, 25 (2), 86–98 .

Messac, A. , & Mattson, C. A. (2004). Normal constraint method with guarantee of even representation of complete Pareto frontier. AIAA Journal, 42 (10), 2101–2111 .

Motta, R. d. S. , Afonso, S. M. B. , & Lyra, P. R. M. (2012). A modified NBI and NC method for the solution of N -multiobjective optimization problems. Structural and Multidisciplinary Optimization, 46 (2), 239–259 .

Mueller-Gritschneder, D. , Graeb, H. , & Schlichtmann, U. (2009). A successive ap- proach to compute the bounded Pareto front of practical multiobjective opti- mization problems. SIAM Journal on Optimization, 20 (2), 915–934 .

Naccache, P. (1978). Connectedness of the set of nondominated outcomes in mul- ticriteria optimization. Journal of Optimization Theory and Applications, 25 (3), 459–467 .

Rockafellar, R. T. (1970). Convex analysis . Princeton NJ: Princeton University Press .

Ruzika, S. , & Wiecek, M. M. (2005). Survey paper approximation methods in mul- tiobjective programming. Journal of optimization theory and applications, 126 (3), 473–501 .

Sanchis, J. , Martínez, M. , Blasco, X. , & Salcedo, J. V. (2008). A new perspective on multiobjective optimization by enhanced normalized normal constraint method.

Structural and Multidisciplinary Optimization, 36 (5), 537–546 .

Sayin, S. (20 0 0). Measuring the quality of discrete representations of efficient sets in multiple objective mathematical programming. Mathematical Programming, 87 (3), 543–560 .

Sayın, S. (2003). A procedure to find discrete representations of the efficient set with specified coverage errors. Operations Research, 51 (3), 427–436 .

Shao, L. , & Ehrgott, M. (2007). Finding representative nondominated points in multiobjective linear programming. In Proceedings of the 1st IEEE sympo- sium on computational intelligence in multicriteria decision making, MCDM 2007 (pp. 245–252). IEEE .

Shao, L. , & Ehrgott, M. (2008). Approximately solving multiobjective linear pro- grammes in objective space and an application in radiotherapy treatment plan- ning. Mathematical Methods of Operations Research, 68 (2), 257–276 .

Vaz, D. , Paquete, L. , Fonseca, C. , Klamroth, K. , & Stiglmayr, M. (2015). Representa- tion of the non-dominated set in biobjective discrete optimization.. Computers & Operations Research, 63 (1), 172–186 .

White, D. (1991). The maximal dispersion problem and the “first point outside the neighbourhood heuristic”. Computers & Operations Research, 18 , 43–50 .

Yu, P. L. , & Zeleny, M. (1975). The set of all nondominated solutions in linear cases and a multicriteria simplex method. Journal of Mathematical Analysis and Appli- cations, 49 , 430–468 .