Univariate global optimization with multiextremal non-differentiable constraints without penalty functions

Uni v ariate global optimization with multiextremal non-dif ferentiable constraint s without penalty functions ∗ Y arosla v D. Sergey ev D.E.I.S. – Universit ` a della C alabr ia, 87036 Rende (CS) – Italy and University of Nizhni Novgoro d, Gagarin A v ., 23, Nizhni Novgorod – R ussia (yaro@si.dei s.unical.it ) Abstract This paper proposes a ne w algorithm for solving constrained global optimiza- tion problems where both the o bjectiv e function a nd constraints are one-dimensional non-dif ferentiable multiextremal Lipsch itz functions. Multiextremal constraints can lead to comp lex feasible region s being collections of isolated points and inter- v als h aving p ositiv e len gths. The case is considered where the o rder the constraints are ev aluated is fi xed by the nature of the problem and a constraint i is defined only ov er the set where the constraint i − 1 is satisfi ed. The ob jective function is defined only over the set where all the c onstraints are satisfied. In con trast to traditional ap- proaches, the new algorithm does not use any additional parameter or v ariable. All the constraints are not ev aluated during ev ery iteration of the algorithm providing a significant acceleration of the search. T he new algorithm either finds lower and upper bounds for the global optimum or establishes that the problem is infeasible. Con ver gence prop erties and numerical exp eriments sho wing a nice performa nce of the ne w method in comparison with the penalty approach are giv en. K ey W ords : Global optimization, multiextremal constraints, Lipschitz functions, continuo us index functions. ∗ This re search was support ed by the follo wing grants: FIRB RB NE01WBBB, FIRB RB A U01JYPN, and RFBR 04-01-00455-a. The autho r thanks Prof. D. Grimal di for proposi ng the application discussed in the paper . 1 1 Introd uction In last decades u niv ariate glob al optimization problems were stu died intensi vely (see [7, 15, 18, 19, 22, 2 8, 3 1, 37, 43]) becau se ther e exists a large n umber of real-life applications where it is necessary to solve such problem s (see [6, 15, 27, 30, 33, 37, 40]). On the other hand , it is imp ortant to stud y these pr oblems because mathe matical approa ches d ev elope d to solve them c an b e gen eralized to the mu ltidimensional case by nu merous schemes ( see, fo r example, on e-point based, d iagonal, simplicial, space- filling curves, and other popular approache s in [14, 16, 17, 23, 26, 29, 37]). Electrotechn ics and electronics are amon g the fields where one-d imensional global optimization method s can be used successfully (see [8, 9, 10, 11, 2 4, 33, 4 0]). Let us consider, for example, the following so-called ‘mask p roblem’ f or tran smitters. W e have a transmitter (for instance, that of GSM cellular phones) that in a frequency inter - val [ a , b ] has an amplitude A ( x ) that should be inside the mask defined by functions l ( x ) and u ( x ) , i.e., it shou ld b e l ( x ) ≤ A ( x ) ≤ u ( x ) . The mask is de fined by international rules agr eed to av oid interf erence app earing when am plitude is too hig h for a given frequen cy and by properties of electron ic componen ts used to construct the transmitter . Then, it is necessary to find a frequency x ∗ ∈ [ a , b ] such th at the power , p ( x ) , of the transmitted signal is maximal. It ha ppens of ten in engine ering optimizatio n pro blems (see [ 37, 4 0]) that if a con - strained is not satisfied at a p oint then many other constraints and the objecti ve function are not defined at this point. This situation holds in o ur mask problem because if f or a frequ ency ξ if hap pens that A ( ξ ) > u ( ξ ) or A ( ξ ) < l ( ξ ) then there is no tr ansmission and the function p ( ξ ) is not defined. Sin ce the amplitude can touch the mask both from its in ternal and its externa l par ts, isolated points in the adm issible region of p ( x ) can take place. If the maximal power p ( x ∗ ) is attained at an isolated point x ∗ , then this point should be discarded from consideration because it cannot be realized in practice. Thus, the solution is acceptable only if it belon gs to a finite interval of a certain length. This p roblem can be refor mulated in th e f ollowing genera l framew ork o f global optimization prob lems considered in this paper . It is necessary to find the g lobal mini- mizers and the global minimum of a function f ( x ) sub ject to constraints g j ( x ) ≤ 0 , 1 ≤ i ≤ m , over an interval [ a , b ] . The objective fu nction f ( x ) and constraints g j ( x ) , 1 ≤ i < m , are multiextremal no n-differentiable ‘black -box’ Lipsch itz functions with a pr iori known Lipschitz constants (to unify the description proce ss the designation g m + 1 ( x ) , f ( x ) is u sed hereinaf ter). V ery o ften in real-life app lications the orde r th e constraints are evaluated is fixed by the n ature of the problem and not all the con straints are de- fined over the who le search r egion [ a , b ] . The worst case is con sidered h ere, i.e., a constraint g j + 1 ( x ) is defined only over subregions where g j ( x ) ≤ 0. This means that if a con straint is not satisfied at a p oint, the re st of co nstraints and the objective function are not defined at that point. Th e sets Q j , 1 ≤ j ≤ m + 1 , can b e so defined as follows Q 1 = [ a , b ] , Q j + 1 = { x ∈ Q j : g j ( x ) ≤ 0 } , 1 ≤ j ≤ m , (1) Q 1 ⊇ Q 2 ⊇ . . . ⊇ Q m ⊇ Q m + 1 . Since the constraints are m ultiextremal, the admissible region Q m + 1 and regions Q j , 1 ≤ j ≤ m , can be collection s of intervals h aving positive lengths and isolated points. Particularly , is olated points appear when one of the constraints touches zero, for exam- ple, if g j ( x ) is the square of some functio n, t hen g j ( x ) ≤ 0 only when g j ( x ) = 0. T o be implementab le in practice, optimal solutions should h av e a f easible neighborh ood of positive leng th thus, an addition al constraint is included in the model: a point x ∗ should 2 belong to an admissible in terval ha ving length eq ual to or g reater than δ > 0 – a value supplied b y the fina l user . The set of a ll such in tervals is d esignated a s Q δ (of co urse, Q δ ⊆ Q m + 1 ). Eventually foun d isolated points and feasible subregion s having len gth less than δ should be exclud ed fro m co nsideration. If the case of infeasible prob lem Q δ = / 0 holds, it should be also determined. W e can n ow state the prob lem for mally . Find the global minim izers x ∗ and the correspo nding v alue f ∗ such that f ∗ = f ( x ∗ ) = min { f ( x ) : x ∈ Q δ } , (2) where the obje cti ve function f ( x ) an d con straints g j ( x ) , 1 ≤ i < m , are multiextremal function s s atisfying the Lip schitz condition in the form | g j ( x ′ ) − g j ( x ′′ ) |≤ L j | x ′ − x ′′ | , x ′ , x ′′ ∈ Q j , 1 ≤ j ≤ m + 1 , (3) and the constants 0 < L j < ∞ , 1 ≤ j ≤ m + 1 , (4) are known (this sup position is classical in global optimiz ation (see [15, 17, 28])). Meth- ods working on the basis of this assumption are called ‘exact’ in literature, metho ds estimating these values are ‘practical’. On the one hand, the exact meth ods ser ve as a basis for studying theoretical properties of practical ones and are used as a unit of measure of the speed of practical metho ds. On the oth er h and, in certain cases, wh en additional inform ation about the objec ti ve function and constrain ts is av ailable, they can be applied directly . An example of such a p roblem is shown in Fig. 1. It ha s two non-d ifferentiable multiextremal constraints g 1 ( x ) and g 2 ( x ) . The corresp onding sets Q 1 = [ a , b ] , Q 2 , and Q 3 are sho wn. Th e point c belongs to the sets Q 1 , Q 2 , and Q 3 but c / ∈ Q δ . Th e set Q δ is shown by the grey color . It can be seen from F ig. 1 that the s ets Q 2 , Q 3 , and Q δ consist of disjoint subregions and Q 2 , Q 3 contain also an isolated point. It is not e asy to find a trad itional algorithm f or solving the problem ( 2)–(4). For example, the penalty appr oach requires that f ( x ) and g i ( x ) , 1 ≤ i ≤ m , are defined over the whole search interv al [ a , b ] . I t seems that missi ng v alues can be simply filled in with either a big nu mber o r the functio n value at the near est f easible poin t. Unfor tunately , in the context of Lipschitz a lgorithms, incorpor ating such id eas can lead to infinitely high Lipschitz constants, causing degeneration of the methods and non-applicab ility o f the penalty approach . A promising appr oach called the index sc heme has be en propo sed in [38] ( see also [39, 40]) in combination with stochastic Bayesian algorithms. An important advantage of the index scheme is that it does not introduce additional variables and/or parameters by op position to classical appro aches in [ 2, 3, 16, 17, 2 5]. It has been re cently shown in [35] that the ind ex sch eme can be also successfu lly used in combin ation with the Branch-an d-Bound appr oach. Unfortun ately , this sch eme ca n not be applied directly for solving the prob lem (2)–(4) because it has good c onv ergence p roper ties when all the sets Q j , 1 ≤ j ≤ m + 1 , have no isolated points – requir ement hardly verified in practice without some additional inform ation about t he pro blem. Thus, isolated points g iv e seriou s pro blems when o ne has only Lip schitz inf orma- tion. First, because it is not po ssible to say a p riori whether th e feasible r egion has isolated poin ts or no t (f or example, the method fr om [35] converges o nly to g lobal minimizers if it is en sured ab sence of isolated points). Second, in Lip schitz glob al optimization isolated p oints ca n lea d to two problem s: accumulatio n of trial points in 3 Figure 1: A n example of the problem (2)–(4). their n eighbo rhood (this happens even if there e xists a sub -region wher e a constraint does not touch zero but is only close to zero ) and increasing the estimates of Lipschitz constants to in finity . This fact means th at the search region will be covered by a uni- form mesh of trials) if the Lipschitz constant is estimated or the method simply will not work if th e L ipschitz constant is given – our case – because local ad aptively obtained informa tion will contradict the g i ven one. Therefo re, traditio nal Lipschitz methods cannot be used in the presence of isolated points and, since their absence can be hardly determined in practice, d ev elopm ents of metho ds that are able to work ind ependen tly of the presence or absence of them becomes very important. In this paper, su ch a method is proposed . It ev olves th e idea of separate consid- eration o f each con straint intro duced in [38] in a ne w way and redu ces th e o riginal constrained pr oblem to a n ew continuou s problem. The meth od from [3 5] is used a s a basis for co nstruction of th e new sche me. Instead of d iscontinuo us sup port functio ns propo sed in [3 5], new continuous f unctions are built. These new stru ctures are very importan t bec ause b y using them it becomes possible to app ly numero us to ols devel- oped in Lipschitz un constrained optimizatio n to a very ge neral class of co nstrained problem s. It is also nece ssary to emp hasize th at th e n ew appro ach do es n ot introduc e additional variables an d/or parameters during this passage from initial discontinuou s constrained partially defined problem to the continuo us uncon strained one. T o conclude this introductio n it is necessary to empha size once again t hat the prob - lem of multi-dimen sional extensions of one-d imensional Lipschitz global optimization methods to many dim ensions is a n on-trivial seriou s pr oblem (P . Hansen and B. Jau - mard write in their su rvey o n Lipschitz o ptimization [15] published in the Han dboo k of Globa l Optimiza tion: ‘Large problems (with 1 0 variables or mo re) app ear to b e of- 4 ten in tractable, at least if high precision is r equired’ ) a nd is beyond th e scop e of this paper dedica ted to the un iv ariate algorithms and un i variate applications. Howe ver , the approa ch p roposed here is very promising from this point of view . I n the futur e, a num- ber of various multi-dimensional exten sions (starting from th e adaptiv e d iagonal and space-filling curves approache s (see [34, 40])) of the alg orithm presented in this pape r will be studied. The rest of the paper is organized as follo ws. The new method is described in Section 2. Section 3 contains comp utational results and a brief conclusion. 2 Continuous index functions and the new algorithm The ind ex scheme (see [ 38, 3 9, 4 0]) co nsiders con straints o ne at a time a t every po int where it h as been decided to calc ulate f ( x ) d etermining the index ν = ν ( x ) , 1 ≤ ν ≤ m + 1 , b y the following condition s g j ( x ) ≤ 0 , 1 ≤ j ≤ ν − 1 , g ν ( x ) > 0 , (5) where fo r ν = m + 1 the last inequality is omitted. The term trial u sed herein after means d eterminin g th e ind ex ν ( x ) at a po int x by evaluation g i ( x ) , 1 ≤ i ≤ ν ( x ) . The index ν ( x ) and the value g ν ( x ) ( x ) are called r esults o f the trial . The disco ntinuou s in dex function J ( x ) , x ∈ [ a , b ] , c an be written for the prob lem (2)–(4) following [38] J ( x ) = g ν ( x ) ( x ) −  0 , ν ( x ) < m + 1 , f ∗ , ν ( x ) = m + 1 , (6) where the value f ∗ is the unknown solution to this prob lem. Let us start our th eoretical co nsideration by noticing that the g lobal minim izer of the original constrain ed problem (2)–(4) in the case Q δ 6 = / 0 coincides with the solution to the following unco nstrained discontinuous problem J ( x ∗ ) = min { J ( x ) : x ∈ Q δ } , (7) where Q δ = [ a , b ] \ { Q m + 1 \ Q δ } . (8) Suppose now that trials have been executed in a way at some points a = x 0 < x 1 < . . . < x i < . . . < x k = b (9) and ν i = ν ( x i ) , 0 ≤ i ≤ k , are their starting ind exes . No te that the notion of index is different with respect to [ 38, 39, 4 0] whe re the in dex is calculate d once and then u sed in the course of optimization. In this pap er , fo rmula (5) defines the starting v alue for the index that can then be changed during the w ork of the algor ithm. The points from (9) fo rm the list (called her einafter History List H ( k ) ) of intervals [ l i , r i ] , 1 ≤ i ≤ k , wher e l i < r i , 1 ≤ i ≤ k , r i = l i + 1 , 1 ≤ i < k . The record x ∈ H ( k ) means th at the point x = l i or x = r i for an interval i from H ( k ) . Every element i , 1 ≤ i ≤ k , of the list contains the following infor mation: [ l i , r i ] , ν ( l i ) , ν ( r i ) , g ν ( l i ) ( l i ) , g ν ( r i ) ( r i ) . (10) 5 The s econ d list, W ( k ) , called W orking List is built durin g the work of the method to be introdu ced by exclud ing fro m H ( k ) inter vals where glo bal m inimizers of the prob lem (2)–(4) can not be lo cated (initially it is s tated W ( k ) = H ( k ) ). In contrast to H ( k ) where the information (10) is calculated once and then is k ept during the search, indexes ν ( l i ) and ν ( r i ) in W ( k ) can be changed in the course of optimization. In order to pass from the problem (2)–(4) to the problem (7) it is n ecessary to estimate the value f ∗ from (2) and the set Q δ . Using th e results of trials at the po ints from the row (9) the value Z ∗ k = min { g m + 1 ( x ) : ν ( x ) = m + 1 , x ∈ W ( k ) } . (11) estimating f ∗ can be calculated if there exist points x with the index ν ( x ) = m + 1 . Th is value allo ws us to define the fu nction J k ( x ) , x ∈ [ a , b ] , by replacin g the unknown v alue f ∗ in (6) by Z ∗ k : J k ( x ) = g ν ( x ) ( x ) −  0 , ν ( x ) < m + 1 , Z ∗ k , ν ( x ) = m + 1 . (12) The follo wing Lemm a establishes some usefu l pr operties of the fu nctions J ( x ) and J k ( x ) . Lemma 1 The following assertions hold for the function s J ( x ) and J k ( x ) : i. for all po ints x having indexes ν ( x ) < m + 1 , it follows J k ( x ) = J ( x ) > 0 ; ii. J k ( x ) ≤ 0 , x ∈ { x : g m + 1 ( x ) ≤ Z ∗ k } . (13) iii. if ν ( x ) = m + 1 and Z ∗ k ≥ f ∗ then J k ( x ) ≤ J ( x ) . (14) Proof . T ruth o f assertions i – iii follows from definitions of the fu nctions J ( x ) and J k ( x ) . Particularly , it follo ws fr om Lemma 1 that (14) hold s if the trial poin t cor respond ing to Z ∗ k belongs to Q δ . Th e estima te (1 4) is no t tru e if Z ∗ k ≤ f ∗ , the situation which can occur only if x ∗ k belongs to the set Q m + 1 \ Q δ where x ∗ k = arg min { g m + 1 ( x ) : ν ( x ) = m + 1 , x ∈ W ( k ) } . (15) Let us intro duce the following con tinuous ind ex functio n C ( x ) , x ∈ Q δ , and stud y its proper ties. C ( x ) = max y ∈ Q δ { J ( x ) , J ( y ) − K ν ( y ) | x − y | ) } , (16 ) where K ν ( y ) such that L ν ( y ) < K ν ( y ) < ∞ is an overestimate of the Lipschitz constan t correspo nding to the fun ction g ν ( y ) ( y ) a nd J ( x ) is the discontinuou s ind ex fun ction from (6). As a n illustra tion, the function C ( x ) corr espondin g to the p roblem pr esented in Fig. 1 is shown in Fig. 2. The parts of the functio n C ( x ) correspo nding to x and y such that J ( x ) < J ( y ) − K ν ( y ) | x − y | are shown by the thin line. 6 Figure 2: T he function C ( x ) correspo nding to the prob lem presented in Fig. 1. If Q δ 6 = / 0 , th e glo bal minim izers of th e or iginal con strained prob lem (2)–(4) c oin- cide with the solutions of the following continuous prob lem C ( x ∗ ) = min { C ( x ) : x ∈ Q δ } . (17) In the case Q δ = / 0 the set Q δ = [ a , b ] \ Q m + 1 and we ha ve ν ( x ) < m + 1 , x ∈ Q δ . Th us, due to Lemma 1, it follows C ( x ) > 0 , x ∈ Q δ . (18) Similarly to d efinition of the f unction J k ( x ) , the value Z ∗ k is used to d efine the function s C k ( x ) , x ∈ [ a , b ] , as follows. C k ( x ) = max y ∈ [ a , b ] { J k ( x ) , J k ( y ) − K ν ( y ) | x − y | ) } . (19) Lemma 2 The following assertions hold for the function s C ( x ) and C k ( x ) : i. ine qualities C ( x ) ≥ J ( x ) , C k ( x ) ≥ J k ( x ) hold over the set Q δ ; ii. if ν ( x ) < m + 1 then C k ( x ) > 0 ; iii. if ν ( x ) = m + 1 and x ∈ W ( k ) then C k ( x ) ≥ 0 ; iv . if x ∈ { x : g m + 1 ( x ) ≤ Z ∗ k } then C k ( x ) ≤ 0 . v . if x ∗ k ∈ Q δ then C k ( x ) ≤ C ( x ) , x ∈ Q δ . Proof . The truth of the assertions follo ws from Lemma 1 and formu lae (11),(15), (16), and (19). It f ollows from L emma 2 that if x ∗ k ∈ Q δ , the global m inimizers can not b e loc ated in zo nes where C k ( x ) > 0 , x ∈ Q δ . Over every in terval [ l i , r i ] we ar e in terested in sub - regions ha ving the index greater or equal to ν i = max { ν ( l i ) , ν ( r i ) } , 7 because, due to construction of the functio n C k ( x ) , only these subregions can probab ly contain a global minimizer . It can be shown (see [35]) that [ l i , r i ] ∩ Q ν i ⊆    [ y − i , y + i ] , ν ( l i ) = ν ( r i ) , [ y − i , r i ] , ν ( l i ) < ν ( r i ) , [ l i , y + i ] , ν ( l i ) > ν ( r i ) , (20) where y − i = l i + z ( l i ) / K ν ( l i ) , (21) y + i = r i − z ( r i ) / K ν ( r i ) , (22) and z ( x ) = J k ( x ) . Let us call any v alue R i , 1 ≤ i ≤ k , characteristic of the interv al [ l i , r i ] if the following inequ ality is true min { C k ( x ) : x ∈ [ l i , r i ] , ν ( x ) = ν i } ≥ R i . (23) It follows from assertion i of Lemma 2 and [35] that (23) is fulfilled fo r R i = ˇ R i where ˇ R i = ˇ R ( l i , r i ) =    0 . 5 ( z ( l i ) + z ( r i ) − K ν ( r i ) ( r i − l i )) , ν ( l i ) = ν ( r i ) , z ( r i ) − K ν ( r i ) ( r i − y − i ) , ν ( l i ) < ν ( r i ) , z ( l i ) − K ν ( l i ) ( y + i − l i ) , ν ( l i ) > ν ( r i ) . (24) The characteristic ˇ R i from (24) depends only on the v alues of the func tion J k ( x ) e valu- ated at the points l i and r i . It does not use any i nfo rmation from other intervals belong- ing to the working list W ( k ) . W e are ready now to introduce the Algorithm workin g with Continuous Ind ex Func- tions (A CIF). It either solves the problem ( 17) or determ ines that the case (18) takes place. Th e ACIF work s by calculating ch aracteristics R i initially using ( 24) an d then improving them durin g th e search by construc ting the function C k ( x ) . On the o ne hand, th e method tries to find a good e stimate Z ∗ k . On th e other han d, it sear ches a nd eliminates from W ( k ) intervals that cannot contain x ∗ using the fact following from Lemma 2 and (23) that if x ∗ k ∈ Q δ , an interval [ l j , r j ] having a characteristic R j > 0 can be elim inated from co nsideration. Th e constrain t introdu cing the par ameter δ h elps to exclude more intervals . Let us take a gen eric interval [ l t , r t ] , 1 ≤ t ≤ q ( k + 1 ) , fr om the working list and calculate its ch aracteristic R ( l t , r t ) . W e will also show how the fu nction C k ( x ) allows us to improve character istics of intervals adjacent to [ l t , r t ] . Initially characteristic for the inter val [ l t , r t ] is c alculated as R ( l t , r t ) = ˇ R ( l t , r t ) . If R ( l t , r t ) ≤ 0 or ν ( l t ) = ν ( r t ) , then th e characte ristic R ( l t , r t ) has b een compute d. If R ( l t , r t ) > 0 and ν ( l t ) < ν ( r t ) go to the operation Bac kwar d motion. Otherwise e xecute the operation Onwar d motion. Backwar d motion. Exclude from W ( k + 1 ) all the inter vals i such that z ( r t ) − K ν ( r t ) ( r t − l i )) > 0 , 1 ≤ j + 1 ≤ i ≤ t − 1 , (25) where the interval j violates (25). Calcu late the value R − j =    0 . 5 ( z ( l j ) + z − ( r j ) − K ν ( r t ) ( r j − l j )) , ν ( l j ) = ν ( r t ) z − ( r j ) − K ν ( r j ) ( r j − l j − z ( l j ) / K ν ( l j ) ) , ν ( l j ) < ν ( r t ) z ( l j ) − K ν ( l j ) ( r j − l j − z − ( r j ) / K ν ( r t ) ) , ν ( l j ) > ν ( r t ) (26) where z − ( r j ) = z ( r t ) − K ν ( r t ) ( r t − r j )) . (27) 8 Figure 3: I mproving characteristics by the operation Bac kwar d motion . If R − j < R j , set in the working list W ( k + 1 ) z ( r j ) = z − ( r j ) , ν ( r j ) = ν ( r t ) , R j = R − j , maintaining in the history list H ( k + 1 ) the original values of g ν ( r j ) ( r j ) an d ν ( r j ) . Cal- culate the number q ( k + 1 ) of the intervals in W ( k + 1 ) . An illustration to the operation Bac kwar d motion is giv en in Fig. 3. Th ree interv als are presented in Fig. 3: [ l i − 2 , r i − 2 ] = [ p , q ] , [ l i − 1 , r i − 1 ] = [ q , h ] , [ l i , r i ] = [ h , d ] . Suppose that the functio n C k ( x ) has been e valuated at the points p , q , h , and d and ν ( p ) = ν ( q ) = ν ( h ) < ν ( d ) . Characteristic ˇ R i − 2 of the interv al [ p , q ] is negative an d the b old lin e sho ws the zone where th e g lobal min imizer co uld be pr obably fo und. The sam e situation h olds fo r the interval [ q , h ] . Since the characteristic ˇ R i of the inter val [ h , d ] is p ositiv e an d ν ( h ) < ν ( d ) , th e operatio n Bac kward motion starts to work . It can be seen fro m Fig. 3 that the new characteristics R i − 2 and R i − 1 calculated using informa tion obtained at the point d are positi ve and , theref ore, the intervals [ p , q ] and [ q , h ] c annot contain glo bal minimizers and can be so excluded from the w ork ing list. Onwar d motion. E xclude from W ( k + 1 ) all th e intervals i such that z ( l t ) − K ν ( l t ) ( r i − l t )) > 0 , t + 1 ≤ i ≤ j − 1 ≤ q ( k ) , (28) where the interval j violates (28). Calcu late the value R + j =    0 . 5 ( z + ( l j ) + z ( r j ) − K ν ( r j ) ( r j − l j )) , ν ( l t ) = ν ( r j ) z ( r j ) − K ν ( r j ) ( r j − l j − z + ( l j ) / K ν ( l t ) ) , ν ( l t ) < ν ( r j ) z + ( l j ) − K ν ( l t ) ( r j − l j − z ( r j ) / K ν ( r j ) ) , ν ( l t ) > ν ( r j ) (29) 9 where z + ( l j ) = z ( l t ) − K ν ( l t ) ( l j − l t )) . (30) If R + j < R j , set in the working list W ( k + 1 ) z ( l j ) = z + ( l j ) , ν ( l j ) = ν ( l t ) , R j = R + j , maintaining in the history list H ( k + 1 ) the o riginal values of g ν ( l j ) ( l j ) and ν ( l j ) . Cal- culate the number q ( k + 1 ) of the intervals in W ( k + 1 ) . In or der to describe the method we need so me d efinitions and initial settings. I t is supposed that: – the search accuracy 0 < ε ≤ δ has been chosen , where δ is from (2); – two initial trials ha ve been executed at the points x 0 = a and x 1 = b ; – I t h as b een assigned W ( 1 ) = H ( 1 ) = [ x 0 , x 1 ] and th e n umber t of the interval to be subdivided at the next iteration has been set to t = 1; – the values Z ∗ k and M k = max { ν ( x i ) : 0 ≤ i ≤ k } (31) have been calculated for k = 1; – the set V δ containing the points x ∈ H ( k ) such that x ∈ Q m + 1 \ Q δ has be en set to V δ = / 0 . Suppose now that k , k ≥ 1 , iteration s ha ve been m ade by the A CIF , the list H ( k ) contains k intervals, W ( k ) contains q ( k ) in tervals fo r which ch aracteristics have b een ev aluated, and an interval [ l t , r t ] for subdivision has been found. T he choice of the ne xt interval to be subd ivided is made as follows. Step 1. (Sub division and the new trial.) Update W ( k + 1 ) and H ( k + 1 ) by substituting the interval [ l t , r t ] in W ( k ) and H ( k ) by the n ew in tervals [ l t , x k + 1 ] , [ x k + 1 , r t ] where x k + 1 =    0 . 5 ( y − t + y + t ) , ν ( l t ) = ν ( r t ) 0 . 5 ( y − t + r t ) , ν ( l t ) < ν ( r t ) 0 . 5 ( l t + y + t ) , ν ( l t ) > ν ( r t ) (32) Execute the ( k + 1 ) -th tr ial at th e po int x k + 1 and, as the resu lt, ob tain the values ν ( x k + 1 ) and g ν ( x k + 1 ) ( x k + 1 ) . Recalculate M k + 1 . Step 2. (Ca lculation of the estimate Z ∗ k + 1 and characteristics.) Associa te with the point x k + 1 the value z k + 1 = J k + 1 ( x k + 1 ) and recalculate the estimate Z ∗ k + 1 if ν ( x k + 1 ) = m + 1 . If Z ∗ k + 1 < Z ∗ k then for all points x ∈ W ( k + 1 ) , x 6 = x k + 1 , such that ν ( x ) = m + 1 set z ( x ) = z ( x ) + Z ∗ k − Z ∗ k + 1 and recalculate c haracterisitcs of the intervals in W ( k + 1 ) . Other wise calculate char acteristics only f or the intervals [ l t , x k + 1 ] and [ x k + 1 , r t ] . Go to Step 3. Step 3. (F ind ing an interval for the ne xt subdivision .) If W ( k + 1 ) = / 0 , then Stop (the feasible region is empty). Otherwise, find in the working list W ( k + 1 ) an interval [ l t , r t ] such that t = min { arg min { R i : 1 ≤ i ≤ q ( k + 1 ) }} (33) and go to Step 4. 10 Step 4. (V erifying ap purtenan ce to the set Q δ .) If the inter val to be subdi- vided can belong to the set Q δ then go to Step 5. Other wise exclude all foun d intervals that are out of Q δ from the workin g list and include the p oints form ing these intervals and ha ving the index m + 1 in the set V δ . If the point x ∗ k + 1 belongs to one of the excluded interv als then go to Step 6 otherwise go to Step 3. Step 5. (V erifying accuracy .) If the inequ ality r t − l t > ε (34) holds, then go to Step 1. In the opposite case, Stop (the req uired accura cy has been reached). Step 6 . (Restarting.) Recalculate the e stimate Z ∗ k + 1 without u sage of the points included in V δ . Form the new set W ( k + 1 ) inclu ding in it all the inter- vals fr om H ( k + 1 ) that do not co ntain points fr om V δ and intervals co ntaining points x ∈ V δ such that z ( x ) > Z ∗ k + 1 . F or all intervals in W ( k + 1 ) rec alculate characteristics R i applying b ackward motio n for a ll intervals i ha ving R i > 0 if ν ( l t ) < ν ( r t ) and on ward motion if ν ( l t ) > ν ( r t ) . In th e latter case, ch aracteris- tics of the intervals satisfy ing (28) are not calculated. Exclude from W ( k + 1 ) all the intervals having positi ve characteristics. Then go to Step 4. Step 4 executes a n important op eration – verify ing ap purtenan ce to the set Q δ . T o do this we ch eck whether the inter val [ l t , r t ] chosen for subdivision can contain a feasible interval having a length greater than δ . Four cases should be considered. i. (Case ν ( l t ) < m + 1 , ν ( r t ) < m + 1 .) I f y + t − y − t = r t − l t − z ( r t ) / K ν ( r t ) − z ( l t ) / K ν ( l t ) < δ (35) then [ l t , r t ] / ∈ Q δ because over [ l t , r t ] only th e in terval [ y − t , y + t ] can possibly co n- tain a global minimizer but its length is less than δ . ii. (Case ν ( l t ) = m + 1 , ν ( r t ) < m + 1 .) An alogou sly , if r t − l t − z ( r t ) / K ν ( r t ) > δ then the interval [ l t , r t ] can belong to the set Q δ . Otherwise, if in the h istory list H ( k + 1 ) there exists an interval [ l j , r j ] , j < t , such tha t ν ( l j ) < m + 1 and r t − l j − z ( r t ) / K ν ( r t ) − z ( l j ) / K ν ( l j ) < δ ( 36) or ν ( l j ) = m + 1 , j = 1 , and r t − l j − z ( r t ) / K ν ( r t ) < δ (37) then all the intervals [ l j , r j ] , . . . , [ l t , r t ] / ∈ Q δ and the correspo nding points r j , . . . , l t ∈ Q m + 1 \ Q δ . iii. (Case ν ( l t ) < m + 1 , ν ( r t ) = m + 1 .) This case is considere d analogously to the previous one but con firmation of possibility fo r [ l t , r t ] to belong to Q δ is searched among intervals i > t . iv . ( Case ν ( l t ) = ν ( r t ) = m + 1 . ) Th is case is a combination of the cases ii and iii. 11 The in troduc ed procedu re verifies inc lusion [ l t , r t ] ∈ Q δ for all possible combin a- tions of in dexes ν ( l t ) , ν ( r t ) . Of co urse, it is a lso possible to simplify Step 4 and verify only condition (35) – the rule determining during the search the major part of intervals belongin g to Q m + 1 \ Q δ . In this case, after satisfying the stoppin g rule from Step 5, it is necessary to ch eck whether the foun d solutio n x ∗ k + 1 belongs to Q m + 1 \ Q δ and, if necessary , t o reiterate the metho d starting from Step 6. The following situations can, therefo re, hold after fulfillment of the stop ping rule: i. Th e alg orithm has finished its work an d the workin g list is em pty , th en Q δ = / 0 and the set V δ contains the points from Q m + 1 \ Q δ if any . ii. Th e workin g list is n ot empty and it d oes n ot contain intervals [ l p , r p ] such that R p < 0 and max { ν ( l p ) , ν ( r p ) } < m + 1 . (38) In this case it is n ecessary to check lo cally in the neighbor hood of x ∗ k whether x ∗ k ∈ Q δ . If this situatio n h olds, then the g lobal min imum z ∗ of the pro blem (2)–(4) can be bound ed as follows z ∗ ∈ [ R t ( k ) + Z ∗ k , Z ∗ k ] where R t ( k ) is the char acteristic corr espondin g to the interval number t = t ( k ) from (33). In the opposite c ase it is nece ssary to inclu de th e p oint x ∗ k in V δ and to return to Step 6. iii. Th e last case considers the situation where the working list is not empty an d t her e exists an interval [ l p , r p ] such that R p < 0 and ( 38) h olds. Again , it is necessary to ch eck locally in the neighborho od of x ∗ k whether x ∗ k ∈ Q δ . If this analysis shows th at x ∗ k / ∈ Q δ then it is n ecessary to include the p oint x ∗ k in V δ and return to Step 6. Othe rwise, the value Z ∗ k can be taken as an up per boun d of the g lobal minimum z ∗ . A lo wer bound can be calculated easily by taking from the w ork ing list the trial points x i such that ν ( x i ) = m + 1 and co nstructing for f ( x ) the support function of the type [28] using only these points. The g lobal minimum of this support fun ction over the in tervals belong ing to the working list will be a lower bound for z ∗ . Consider no w the infinite trial seq uence { x k } generated by the algorithm ACIF when ε = 0 in the stopping rule (34). W e denote by X ∗ the set of the global minimizers of the problem (2)–(4) and by X ′ the set of limit points of the sequenc e { x k } . The following two theor ems d escribe conv ergence con ditions of the A CIF . Since they can be d erived as a particular case of gener al convergence studies gi ven in [17] (Branch - and-Boun d approach) and [32] (Divide the Best algorithms) their proofs are omitted. Theorem 1 If the pr o blem (2)–(4) is feasible, i.e . Q δ 6 = / 0 , then X ∗ = X ′ . Theorem 2 If the p r o blem (2)–(4) is infeasible then the algorithm ACIF stops in a finite number of iterations. 3 Numerical comparison and conclusion The A CIF has b een numerically compa red to the algorith m (in dicated here inafter as PEN ) proposed b y Pijavskii (see [28 , 15]) combined with a penalty f unction. The 12 T able 1: Nu merical results o btained by the PEN on 10 no n-differentiable and 10 d if- ferentiable proble ms. Problem Non-dif ferentiable Differentiab le Iterations Eva luations Iterations Evaluations 1 247 494 83 166 2 241 482 954 1906 3 797 1594 119 238 4 272 819 1762 5286 5 671 2013 765 2295 6 909 2727 477 1431 7 199 597 917 2751 8 365 1460 821 3284 9 1183 4732 262 1048 10 135 540 2019 8076 A verage 501 . 9 1545 . 8 817 . 9 2648 . 1 PEN has been chosen for compa rison because the method of Pija vskii in l iteratur e (see [14, 1 6, 17, 23, 26, 29, 37]) is used as a kind of the unit of measure of efficiency of the ne w Lipschitz glob al op timization algorithms and it uses in its work the sam e informa tion about th e pr oblem as th e ACIF – the Lipschitz co nstants fo r the o bjective function an d constrain ts. The u sage o f the penalty scheme allows us to emph asize advantages of the index approac h. Since th e PEN in every itera tion ev aluates the objective fu nction f ( x ) and a ll the constraints, twenty feasible test p roblems (ten d ifferentiable an d ten n on-differentiable) introdu ced in [13] h av e been used for testing th e new algorithm. The A CIF has also been app lied to on e d ifferentiable an d one non-differentiable infeasible test problems from [13]. In all the e xper iments there has been ch osen the or iginal ( see [13]) ord er the constrain ts are ev aluated du ring optimizatio n, without d etermining the best f or the A CIF order . In th e PEN, the co nstrained problems were re duced to the unco nstrained on es as follows P ∗ ( x ) = f ( x ) + P max { g 1 ( x ) , g 2 ( x ) , . . . , g N v ( x ) , 0 } (39) and coef ficients P from [ 13] h av e been used. The same accur acy ε = 10 − 4 ( b − a ) (where b and a ar e fr om (2)) and the star ting tr ial po ints a and b have been used in all the experiments for both A CIF and PEN. T able 1 contains numerical results obtained for the PEN. The column “Ev aluation” shows the total num ber of e valuations equal to ( N v + 1 ) × N iter , where N v is the num ber of con straints and N iter is the number o f iteration s for each problem . T ables 2 and 3 present numerical results for the new method for δ = ε and δ = 1 0 ε . The columns in the T ables hav e the following meaning : - th e colu mns N g 1 , N g 2 , an d N g 3 present th e nu mber o f trials where the con straint g i , 1 ≤ i ≤ 3 , was the last e valuated constrain t; - th e column N f shows ho w many tim es the objecti ve function f ( x ) has been ev al- uated; 13 T able 2: Results ob tained by the new algorith m on the non- differentiable problems. Problem δ = ε δ = 10 ε N g 1 N g 2 N g 3 N f Iter . Eva l. N g 1 N g 2 N g 3 N f Iter . Eva l. 1 23 − − 28 51 79 23 − − 28 51 79 2 18 − − 16 34 50 17 − − 16 33 49 3 95 − − 18 113 1 31 80 − − 18 98 116 4 107 14 − 84 205 387 82 11 − 84 177 356 5 153 88 − 24 265 401 114 66 − 24 204 318 6 16 16 − 597 629 1839 16 15 − 597 628 1 837 7 52 18 − 39 109 205 49 14 − 39 102 194 8 28 11 3 21 63 143 28 11 3 21 63 143 9 8 81 49 183 321 1049 8 59 32 183 282 954 10 32 3 17 13 65 141 30 2 17 13 62 137 A verage 53 . 2 33 . 0 23 . 0 102 . 3 185 . 5 442 . 5 44 . 7 25 . 4 17 . 3 102 . 3 170 . 0 418 . 3 T able 3: Results ob tained by the new algorith m on the differentiable problems. Problem δ = ε δ = 10 ε N g 1 N g 2 N g 3 N f Iter . Eva l. N g 1 N g 2 N g 3 N f Iter . Eva l. 1 10 − − 13 23 36 10 − − 13 23 36 2 199 − − 21 220 241 155 − − 21 176 197 3 40 − − 22 62 84 38 − − 22 60 82 4 480 127 − 189 796 1301 212 73 − 189 474 925 5 8 13 − 122 143 400 8 13 − 122 143 400 6 14 55 − 18 87 178 13 34 − 18 65 135 7 36 13 − 241 290 785 35 13 − 241 289 784 8 94 21 5 82 202 479 80 19 5 82 186 461 9 7 35 6 51 99 299 7 32 6 51 96 293 10 36 14 174 1173 1397 5278 35 10 92 1173 1310 5023 A verage 92 . 4 39 . 7 61 . 7 193 . 2 331 . 9 908 . 1 59 . 3 27 . 7 34 . 3 193 . 2 282 . 2 833 . 6 - th e column ”Eval. ” is the total n umber of evaluations of the o bjective f unction and the constraints. T his quantity is equal to: - N g 1 + 2 × N f , for problems with one constraint; - N g 1 + 2 × N g 2 + 3 × N f , for pro blems with tw o constrain ts; - N g 1 + 2 × N g 2 + 3 × N g 3 + 4 × N f , for pro blems wi th th ree constraints. It can be seen from the T ables that in all th e experiments the A CIF significan tly outperf orms the PEN b oth in iterations and ev aluations. The AC IF works faster if th e difference between δ and ε increases. Th is effect is especially notab le for problems where it is necessary to e xecute many iterations out of the feasible re gion (see columns N g 1 , N g 2 , N g 3 for non- differentiable problems 3– 5, 9 an d d ifferentiable problems 2, 4 , 8, 10). Note that the penalty app roach requires an accurate tun ing of the penalty coeffi- cient in contrast to the ACIF that works with out nece ssity to d etermine any add itional parameter . Mor eover , whe n th e pen alty appr oach is used a nd a constraint g ( x ) is de- fined only over a subr egion [ c , d ] o f th e sear ch region [ a , b ] , the prob lem of extending g ( x ) to th e whole r egion [ a , b ] arises. The A CIF does not have this d ifficulty be cause 14 the constrain ts and th e objective function are ev aluated on ly within their r egions of definition. Finally , the penalty app roach is not able to determ ine whether a prob lem is in - feasible. The A CIF with δ = ε has determin ed infeasibility of the non-d ifferentiable problem from [1 3] in 86 iterations consisting of 81 ev aluatio ns of the first constrain t and 5 ev aluations of the first an d second co nstraints (i.e., 91 evaluations in total). Th e infeasibility of th e differentiable prob lem from [13] h as been de termined by the ACIF with δ = ε in 38 iteration s con sisting o f 9 evaluations of the first con straint and 2 9 ev aluations of the fir st a nd secon d constra ints (i.e. , 67 ev aluations in total). Naturally , the objective functions have not been ev aluated in both cases. In conclusio n, we illu strate perf ormance of the new m ethod (see Fig . 4) and the PEN (see Fig. 5) on the non-d ifferentiable problem 9 from [13]. min x ∈ [ 0 , 4 ] f ( x ) = 3 − 2 exp  − 1 2  22 5 − x      sin  π  22 5 − x      subject to g 1 ( x ) = 3  exp  −     sin  5 2 sin  11 5 x       + 1 100 x 2 − 1 2  ≤ 0 , g 2 ( x ) =          6  x − 1 2  2 − 1 2 x ≤ 1 2 1 4  x − 5 2  x > 1 2 ≤ 0 , g 3 ( x ) = 4 5 −      sin  24 5 − x      + 6 25 − x 20  ≤ 0 . The pro blem has 3 disjoint f easible subregions shown in Fig. 4 by co ntinuou s bold in- tervals o n the line f ( x ) = 0, the global optimum is located at the p oint x ∗ = 0 . 95019 236 (see Fig. 4). The objective f unction is shown by a solid lin e and the constraints are drawn by dotted/mix-do tted lines. The first line ( from up to do wn) of “+” lo cated und er the graph of the pr oblem 9 in the upper subplot of Fig. 4 represents the points wh ere the first co nstraint has n ot been satisfied (number of iteration s equal to 8). T hus, due to the de cision rules o f the A CIF , the second con straint has not been evaluated at th ese poin ts. The seco nd line of “+” represents the p oints wher e the fir st con straint has been satisfied but the second constrain t has been not (n umber of iterations eq ual to 59). In th ese points bo th constraints have been ev aluated but the objective function has been no t. The third lin e of “+” r epresents th e po ints wh ere b oth th e first an d th e seco nd co nstraints h av e been satisfied b ut the third constraint has been no t (nu mber of iterations equal to 32). T he last line re presents the points wher e all the constrain ts have be en satisfied and, therefore, the objectiv e fun ction has been ev aluated (n umber o f evaluations equal to 183). The total num ber of ev aluation s is equal to 8 + 59 × 2 + 32 × 3 + 183 × 4 = 954. These ev aluations have been executed durin g 8 + 59 + 32 + 183 = 282 iteratio ns. The lower subplot in Fig. 4 shows dynam ics of the s earch . Fig. 5 sh ows the pena lty function co rrespon ding to P = 15 an d dy namics of th e search e xecuted b y th e PEN. The line of “+” loca ted u nder th e gr aph in the upper subplot o f Fig. 5 represents the points wh ere th e fu nction ( 39) h as b een e valuated. The numb er of iterations is equal to 1 183 and the num ber of ev aluations is equal to 1183 × 4 = 4732 . 15 Figure 4: Behaviour of the ne w method on the non-differentiable problem 9 from [13]. 0 0.5 1 1.5 2 2.5 3 3.5 4 −1 0 1 2 3 x 0 0.5 1 1.5 2 2.5 3 3.5 4 50 100 150 200 250 300 x Nr. of Iterations Refer ences [1] Archetti F . a nd F . Scho en (1984 ), A survey on the glo bal optimiz ation problems: general theory and computational approaches, Anna ls of Operations Resear ch , 1 , 87–11 0. [2] Bertsekas D.P . (1996), C onstrained Optimization and Lagrange Multiplier Meth- ods , Athena Scientific, Belmont, MA. [3] Bertsekas D. P . (19 99), Nonlinea r P r ogramming , Second E dition, Athena Scien - tific, Belmont, MA. [4] Bomze I .M., T . Csendes, R. Hor st, and P .M. Pardalos (1 997) Developments in Global Optimization , Kluwer Academic Publishers, Dordrech t. [5] Breiman L. an d A. Cutler (1 993), A d eterministic algorithm for g lobal optimiza- tion, Math. Pr ogramming , 58 , 179– 199. [6] Brooks S.H. (195 8), Discussion of random methods for locating surface maxima, Operation Resear ch , 6 , 244– 251. [7] Calvin J. and A. ˇ Zilinskas (1999), On the con vergence of the P-algorith m for one-dim ensional global optimization of smooth functio ns, JO T A , 102 , 479–4 95. [8] L. G. Casado, I. Garc ´ ıa, and Y . D. Serge yev (2000 ), I nterval bran ch and boun d algorithm for find ing the Fi rst-Zer o-Crossing-Po int in one-dimensional functions, Reliable Computing , 2 , 179–1 91. 16 Figure 5: Behaviour of the PEN on the non-differentiab le problem 9 from [13]. 0 0.5 1 1.5 2 2.5 3 3.5 4 0 5 10 15 20 25 x P*(x) 0 0.5 1 1.5 2 2.5 3 3.5 4 0 200 400 600 800 1000 1200 x Nr. of Iterations [9] L. G. Casado , I. Garc´ ıa, and Y . D. Sergeyev ( 2002 ), I nterval algorithms for findin g the minimal ro ot in a set of multiextremal no n-differentiable o ne-dim ensional function s, SIAM J. on Scientific C ompu ting , 24(2) , 359–37 6. [10] P . Daponte, D. Grim aldi, A. Mo linaro, an d Y . D. Sergeyev (1995), A n a lgorithm for fin ding th e zer o cr o ssing of time sig nals with Lip schitzean deriva tives , Mea- surements, 16 , 37–4 9. [11] P . Dap onte, D. Grimaldi, A. Mo linaro, and Y . D. Sergeyev (1996) , F ast detection of the first zero-cro ssing in a measurem ent s ignal set, Measu r ements , 19 , 29–39. [12] Evtu shenko Y u. G., M.A. Potapov and V .V . K oro tkich (1992), Numerical methods for global op timization, Re cent Advance s in Globa l Op timization , ed. by C.A. Floudas and P .M. Pardalos, Princeton Uni versity Press, Princeton. [13] Famularo D., Sergeye v Y a.D., and P . Pugliese (2002), T e st Problems for Lipschitz Univ ariate Glo bal Optimization with Multiextrem al Constrain ts, Stochastic an d Global Optimization , ed s. G. Dzemyda, V . Saltenis and A. ˇ Zilinskas, Kluwer Academic Publishers, Dordre cht, 93-110. [14] Floud as C.A. and P .M. Pardalo s (19 96), Sta te of the A rt in Glob al Op timization , Kluwer Academic Publishers, Dordrecht. [15] Hansen P . and B. Jaumard (1995), Lipshitz optimization. In: Ho rst, R., and Parda- los, P .M. (Eds.) . Han dboo k of G lobal Optimization , 407-4 93, Klu wer Academic Publishers, Dordrecht. 17 [16] Horst R. an d P .M. Pardalos (1995), Hand book of Global Optimization , Kluwer Academic Publishers, Dordre cht. [17] Horst R. and H. T uy (1996 ), Global Op timization - Deterministic Appr oaches , Springer–V erlag , Berlin, Third edition. [18] Lamar B.W . ( 1999) , A metho d for converting a class o f u niv ariate function s into d.c. functions, J. of Globa l Optimization , 15 , 55–71. [19] Loca telli M. and F . Schoen (1995), An adapti ve stochastic global optimisation al- gorithm fo r one-d imensional functions, An nals of Op erations r esear ch , 58 , 2 63– 278. [20] Loca telli M. an d F . Scho en ( 1999 ), Random Linkage: a family o f accep- tance/rejection algorithms for global optimisation, Math. Pr ogramming , 85 , 379– 396. [21] Lucid i S. (1994 ), On the role of continu ously differentiable exact pen alty f unc- tions in constrained global optimization, J . of Globa l Optimization , 5 , 49–68. [22] MacLag an D., Sturge, T ., and W .P . Baritom pa (19 96), Equ iv alent Methods for Global Optimiz ation, State of the Art in Glo bal Optimization , e ds. C.A. Flou das, P .M. Pardalos, Kluwer Academic Publishers, Dordrecht, 201–21 2. [23] Mladin eo R. (1992 ), Con vergence rates of a global optimizatio n algorithm, Math. Pr ogramming , 54 , 223 –232. [24] Molin aro A., Sergeyev Y a.D . (200 1) An efficient algorithm for th e zero-c rossing detection in digitized measureme nt signal, Measurement , 30(3) , 187–196. [25] Noced al J. an d S.J. Wright ( 1999 ), Numerical Op timization (Springer Series in Operations Research), Springer V erlag. [26] Pardalos P . M. and J.B. Rosen (19 90), Eds., Computationa l Methods in Global Optimization, Annals of Operations Resear ch , 25 . [27] Patwardhan A.A. , M.N. Karim and R. Shah (1 987), Controller tunin g b y a least- squares method , AIChE J. , 33 , 1735–1737 . [28] Pijavskii S.A. (1972), An Alg orithm for Finding the Abso lute Extrem um of a Function, USSR Comput. Math. and Math. Physics , 12 , 57–67 . [29] Pint ´ er J.D. (1996), Glob al Optimization in Action , Klu wer Academic Publish er , Dordrech t. [30] Ralston P .A.S., K.R. W atson, A.A. Patw ardh an and P .B. Deshp ande (19 85), A computer algorithm for optimized control, Indu strial a nd Engin eering C hemistry , Pr o duct Resear ch and Development , 24 , 1132 . [31] Sergeyev Y a. D. (199 8), Global o ne-dimen sional optimization using smoo th aux- iliary function s, Mathematica l Pr ogramming , 81 , 127-146. [32] Sergeyev Y a.D. (199 9), On con vergen ce of ”Di vide the Best” global optimization algorithm s, Optimization , 44 , 303– 325. 18 [33] Sergeyev Y a.D., P . Daponte, D. Grima ldi and A. Molin aro (1 999), T wo me thods for solving optimization problems arising in electronic measurements and electri- cal engineering , SIAM J. Optimization , 10 , 1–21. [34] Sergeyev Y a.D. (20 00), An Ef ficient Strategy for Adaptive Partition of N- Dimensional Intervals in the Frame work of Diagonal Algorith ms, J ournal of Op- timization Theory and Application s , 107 , 145– 168. [35] Sergeyev Y a. D., Famularo D., and P . Pugliese (2 001) , In dex Bran ch-and -Bound Algorithm for Lip schitz Uni variate Global Optimization with Multiextremal Con- straints, J . of Globa l Optimization , 21 , 317–341 . [36] Sergeyev Y a.D. and D.L. Markin (1995), An algorith m fo r solving global o pti- mization problems with nonlin ear constraints, J . of Global Optimization , 7 , 407– 419. [37] Stron gin R.G. (19 78), Numerical Methods on Multiextr emal Pr oblems , Nauka, Moscow , (In R ussian). [38] Stron gin, R.G. (198 4), Numerical methods for multiextremal nonlinear prog ram- ming prob lems with no nconve x constraints. In: Demy anov , V .F ., an d Pallaschke, D. (Eds.) Lecture Notes in Economics a nd Mathematical Systems 255, 278-282. Proceeding s 1984 . Springer-V erlag. IIASA, Laxenburg/Austria. [39] Stron gin R.G. and D.L. Markin (1986 ), Minimizatio n of multiextremal function s with nonconve x constraints, Cybernetics , 22 , 486–4 93. [40] Stron gin R.G. and Y a.D. Sergeye v (2000), Global Optimization with Non-Con vex Constraints: Sequen tial a nd P arallel Alg orithms , Klu wer Acad emic Publisher s, Dordrech t. [41] Sun X.L. and D. Li (1999), V alue-estimation fu nction metho d for constrained global optimization, JO T A , 10 2 , 385–409 . [42] T ¨ orn A. and A. ˇ Zilinskas (19 89), Global Optimization , Sprin ger–V erlag, Lecture Notes in Computer Science, 350 . [43] W ang X. and T .S.Chang (199 6), An improved uni variate glo bal optimization al- gorithm with improved linear bound ing fu nctions, J. of Global Op timization , 8 , 393–4 11. [44] Zhig ljavsk y A.A. ( 1991) , Theory o f Global Random Sear ch , Kluwer Academic Publishers, Dordrecht. 19