Optimització de Models No Lineals Sense Constriccions

EST APNL APLICACIONS DE LA ,, PROGRAMACIO NO LINEAL Optimització de Models No Lineals Sense Constriccions F. Javier Heredia ' ¡¡ ) ~ - UN I VERSIT AT POLITECNICA DE CAT ALUNYA Biblioteca l l ll l ll lllll lllll lllll lllll llllllll ll lll ll l ll ll ll lll llll /111 1400280220 Optimització de Models No Lineals Sense Constriccions Sudden lnfant Death Syndrome (SIDS) F. Javier Heredia Cervera Departament d 'Estadística i lnvestigació Operativa Secció d'lnformatica UPC Curs 96-97 (111( !') ~- F~cutt.· .• ~ .i , ). • •• ;~¿!".'f',.::~f.! i EüimM'i!1t:.n u f;¡tinoteca Í11cJl'X 1 Índex 1 Sudden Infant D eath Syndrome. . . . . . . . . . . . . ............ 3 l. l. Presenta.ció d<>l prnlikma . . . . . . . . . . . . . . . . . . . . . . . . . . 3 i.2. ;._ro,1c 1 . . . . . . . . . . . . . . . . . ... . .......... . ... . ....... . 3 1.3. Problema d 'o ptimi t.zaciú . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4. Mét.ode ele resolnció . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.5. Realit.zació ele l'exPrcici. . . . . . . . . . . . . . . . . . . . . . . . . . 5 1. 5.1. Gcncracicí de rh1dt> s . . . . . . . . . . . . . . . . . . . . . . 5 1.5.2. Comprnvncic> de les n1tines rl'usiwri . . . . . . . . . . . . . . . . . . . . . . G 1.5.3. Procés d'optimit.z<1ci< 0> • • • • • • • • • • • • • • • • • • • • • • • • • • D 1.5.3.1. Prog;rnma Principal . . . . . . . . . . . . . . . . . . D 1.5.3.2. ;._r unt<ttg;e amli l<'s rntiu es NAG MarklG ........... 11 l.G. Informe dce la pri1.ct.ica .......... . 12 2 Docum entació de la rutina E04LBF . ... ......... .... .. .. 13 2.1. Plaiu T<>xt. DPsnipt.iou . . . . . . . . . . . . . . . .... . .. . ....... . ...... 13 2.1.1. Purpose .. .. .. .. ................ ....... . ...... 13 2.1.2. S¡)f'cifirntion . . . . . .... . ......... .. . . .. . . ... . ... . . 13 2.1.3. Param et.ers .. ... .... ... .................. 13 2.1.4. Erro r Iudir·at.ors <tllll \Vaniiug;s . . . .. . ....................... 21 2.2. Exam ple Prng;ram Somce Crnle . . . . . . . . . . . . . . . . . . . . . . 23 2.2.1. .1\Ia iu prog;rnm . . . . . . . . . . . . . . . 23 2.2.2. 01Jjcctive fnn rr i ou awl firsts rlerivativcs . . . . . . . . . . . . . . . 25 2.2.3. SPco nd de ri ,·;1tiv<'s . . . . . . . . . . . . . . . . . . . . . . 2G 2.2.4. Iufonnation ont¡mt . . . . ... . .. . . .. . .. . . . .. . . . .. ......... 2G 3 1 Sudden Infant Death Syndrome. 1.1 Pr ese nta c ió d el problema L'o bj< 'cti u (Llf¡m·sta pr i1ct ica és r1j11.star 11u matf'ix mo del lrn nno nic a du es ser i es d'ob servn cion s. Lc-1 pl'i111cr::1 t .. sfi 1 fon ua da per lf's 111orts diá1·ies cl ,i uf an t.s catalog;;tdf's com a SIDS al ll<tr g; del prriodP de cinc a 11ys 10i9-83 a Ang;lat. crrn i Gales. La sP, 1..~pna cousist.P ix e11 els Yttlors diaris de l a. ternperatura pro rnig en la matcix a z ona al llnrg; dd matPix 1wriode. Ac¡1w st. ajust és una de les par t.s de l' es tndi d<' :'\Im pliy i CamplH'll so lm~ b rebci<Í C'Ut.re les morts ¡wr SIDS i la temperatura < 1 > 1.2 Model Es di spos<1 d · 11u c onjuut d 'o l>S r ·n· aci011s cid nombre de mort.s cliiL ri es pe r SIDS . fl\IC' iwlicm " :' lll ]Wr lll ;.i. i dP tem¡wrat.ures, f;.i, on l'inclex i recorre els anys (i = 1.2.3.4 ,5) i l'íwkx ./ n·rnn <' ds di es (.j = l, .. ., 3G5). Per t.a! ele simplificar 11 0 tiwl n•111 eu compre!;, ¡n<'s< \u cin. d' un n. ny bixest. dins del pcriode d'estucli, co11~·Ü<lert·111t <1xí qne es disposa clt- ~ 3G5 ohs t;)r vacions pf •r a ca da any . Es vol n justm d rn<1teix 1110<1"1 hiin11!mic" [¡, shie de mort .s i ele tcm¡wrat u res 111 ,1 .i(a¡n 1l1t.ci") = ui 11 ¡ ·os( .. •)· j - ht ) +et t i.1 (u~, h:. r;) ==u~ cos( w.. ' · j - b~) + e~ i = l.2 ,3 . 4,5 j E [l. 3G:i) (1.1) (1.2) <I> :'\Imphy . :'II.F.G .. :'II..J. Ciimpl wll Sl!(/rlc11 i11fr111 t de;Llli syndrorne m1d e11v iro111 11c11tnl te111JH- ' l'i - 1t11n ': cu 1 < 111 c il., -sis· 11sii1g \'ital stst i sti c s, .To 11n u 1l of EpidPniiology ancl Comurn11ity He <i ltli. 19 SI . -11. G3 -- il. 1.3 1:><"-> 3 .':es ~-.;_rs - ~.e:-:·a _:e n ~·=:~~:a e~s .S S : : ) : 2) !"'J - :r:.:o- 4 D('pt..EIO / r. l-'C / O¡it ÍIJJÍI zacú) de rwulels lJO li11ea ls sel!Sf' C'Otist riccions n.mh "' = 2ir /3G3. Aqu <"s ts dos wod<'b de¡wnen de tres pari11u et nes, a;, b¡ i e;, a mb valor difrrt>11ts per rnd<1 ª "Y· Aixú ¡iroporrio1rn un t.ot.al de r¡uinze pariu11etres al lh11g ele ! ¡wr iod t> de ciuc <1 ny s. 1.3 Problema d'optimit za ció L 'a jn st. dt ' les seri<,S d t> mort:; i ck temperatura és fara per separat, donan t !loe a clos prob1"mes d'op t. i mitz aciú diferP 11t. s. Ca.dascun d'e lls con si st ira en la minimit. za ció el e la snma de r edidus entre el valor ob serva t t;j (o m;j) i el valor calcu lat t; (j) (o m , (j) ), al quadrn t. Si plantt'.i<'m < H¡u e,; ta fon ció objed . iu per les du es ser ies qu e es vol ajustar s 'olitenen d prob le mes d'opt.imitz ac ió se güents: ( Pm) ( Ptl llllll f( ./') = 11: 11 , 11: 11 ,c: 11 ' í =i,2 ,:1,4,[1 ¡njn , f(r) (1 ¡ ,/1; ,C¡ i = l, 2 , ~1 .·I J• 1 .~ :H;!"¡ . m m 111 2 :)" LL(111;1(a; , h, ,e, )-m ;1) -i =l .i=l 1 r1 :ti¡!) :)" L L(f;j(a:,11 !,c:J- t;¡) 2 _, i=J ,1 = 1 011 el vecto r de varia.bles del prnblem;1 ( P m) é·s: .1 := (<1; 11 !Jt e~ª ('~ l ) i pel probl ema ( P t) i·s: .r =( u; 1i; ( ' ~ e~) 1.4 Métode de resolució ( 1.3) (1.4) ( 1.5 ) (1.G) Els dos prnblenws (P t) i ( P m) s ' lrnm a .n ele resoklre mit.jan\a.nt la ru tina N AG e 041h f que impknwn ta. t\ll mi•to de el e Nt'wt.on moclificat acbp t. a.t . per a la resolució de probl emes d' op t. im itzacicí 110 linea l a.mb o s ense fit.e s a. l es variable s. Dins d' aquest clon1me11t. es po t tro1"1r la clonunen t. aciú complet a de la ru t ina NAG e04l bf La. forma liabit.ual cl 'ús ck lc·s r ntines cl' op t.imit.za.ció de les lli hr eries 11urn e riques usnal s (Ha.rw<-'11, :\'AG, I:\lSL , etc) rnusi., t<>ix en la prognmrnció, pcr pa.rt ele 1 \1 s uari el 1 u11 pn 1 gro:rr w. ¡rr úu :iprd i d · 1111e s · n1 . t úu ~ .~ rl 'u .. ., ·1w ,r i: • Prng ra.nw. ]11'incipnl: té com ;1 funcions primonli<1ls a.) n ear el mo del que es J Suddeu Info .ut Death SyJJdro111e 5 vol optimitzar, b) inicialitzar ds pariunctrcs de la rutina d'opt.imit .z a.ció, e) frr la crida a la rutina el'optimitzacicí i d) presentar els resultats obt.inguts. • Rutine.q d '·11.crnari: proporeioueu a la rutina el 'opt.imit.zació la informa.- ció neces saria per tal ele resol dre el problema. En el ca.s ele la rutina d'optimit . zac ió e04lbfles rutirn•o cl'usuari són tres: * FUNCT : rutina que, a partir del vector de variables x, calcula el valor de la funció object.iu f(x) i del gradieut Y' f(x). * HESS: a partir del vector de variab les x calcula el valor dels elements ele la mat . rin hessiana \' 1f(.r). * J\IOXIT: rutina de control ele·! procés d'optimització. És cridada pt 'r la rutina d'opt.imit.za cic'> e04l/1f rada IPRINT it.Prn.eious, cm IPRINT és tlll deis pmi1.metres d\,ntrada de e04Jbf L'usuari pot , cada vega.- da. que' es c:rida a .'\IOXIT c•mc't.n • informacic'> sobrf' l'P stat del procés d 'o ptiinit.zació a la itcr;1cic'J <:1 c t.11al. La descripcic\ deta.llacb d'aqnestes t.r<'s rutiues es pot t.robar en la docurne11tació inclusa. e11 acptest. clon111H'llt. 1.5 Realització de l'exercici. 1.5.1 Generació de dad es Cada a.lumne disposa cl'nn coujuut. ele clades {m;j , Í;j; ·i = 1, ... , 5, j = 1, ... ,3G5} cliforents. A<¡nest.<>s dacles s'han de generar mnb l 'e xecuc.ió del programa sidsg¡)J'o/J c¡n<" es trnha iustal·lat al clirectori DIR$EID: [DNLC] tkl Cluster VAX de la FID. Si es fa dese!<" <¡nalsevol tlirectori: RUN DIR$EID: [ONLC] SIDSGPRDB {R.et11rn} apareix per pantalla.: 6 :r:. > - -- > - - . 1 rr::, 1. 5.:? ( 5 e ~ · .;. · _:i_.::_a ~::;s - ~70- ... ·-~-- ... =.e:- 6 Dept.F:IO / UPC / Opli1J1itzaciú de rnodels IJO lilleals se11se coustriccions ... ... .......... .... .. .. . . . . . . . . . . . . . . . . . . . . . . . . . : : ONL, Dept . EIO. UPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . OPTIMITZACIO DE MODELS NO LINEALS SENSE CONSTRICCIONS 'SUDDEN INFANT DEATH SYNDROME (SIDS) ' Cre acio del fitxer de dades pels alumnes -Numero d'identificacio? 100 { Ret11rn} > Dades 100 cr ea de s. > Fitxer sidslOO.dat creat . on els pariunetres en iUtlica ( 100 {Ret11rn}) r eprese nt en les dad es introduid es d es del te rminal per l' alumne. En a.quest exe mpl e s' ha denrnnat la crea.c ió de l es dades corresponents al p roJ,J, ~ ma número 100. sidsgpr ob !mura crea t en el directori actual un fit.xer el e uom sids100 .dat mnb clues columnes de 3 G5 x 5 = 1825 elenwnts, la primPrn . cmTtesponent als valor obs er vat s de morts (m;j) i la se gona a. t.emperatmes (f; 1 ): 4 -4 .0 -->m11=4 f¡¡ = -4.0 6 -1. 5 __, 1n 12 = G f12 = - 1.5 4 -1 . 7 ---> n11:1 = 4 11:1 = -1.7 5 7.7 ___,. ' /H r,,:H i4 = 5 t , ,,] (;.j = 7.7 5 7.4 ___,. 1n.5 ,: tt; .I) = 5 f;),: l!i !j = 7.4 1.5.2 Comprovació de les rutines d'usuari Un co p cr ea d es les rnt.iuf's FL-!\-CT i HESS s' ha de pro ce dir a la seva comprovac ió co m a pa s pr ev i 1wctessa ri a l'ú s de la rutina. cl ' optimit. z ació e04Jb f Per tal de r ealit.zar aquesta comprnv aci< '> s 'lrn cl eixat a cli s po sició dei s alum nes el programa. s id stes t. Aqu est prog;rnma. s 'in voca. fent: RUN DIR$EIO: [ONLC] SIDS TE ST {Re t¡¡rn} L'u snari ha. d'lrnvf'r creat. pn'·, ·i:11w• 1it. nn fit.xer de nom sids. var format per se t. ze lini<'s, la prim Pra. con t.enint. <'l nmnhre ele va ria bl es ( 15) i l es quin ze restants cont . enint el v alor 1111m<~ric de ks \':ui: il,Jes ,¡,.¡ pr obl ema , en l 'or dre indic at a la 1 Sudden Iiif;rnt D('ftf Ji SyJJdnH1u~ 7 formulació (1.4) i (1.3) (c\s a. dir , u1• h¡. c1, ... ,1!5 ,li5, e;,). Es necomana, per t al de no enuuascarar possibles enors, i11trod11'ir 1111 p11nt :r amb component.s cmuplet.ameut difereuts. l'u exemple ele fitxer sids. var seria : 15 1. 1 2.1 3.1 1. 2 2.2 3.2 1. 3 2.3 3.3 1 .4 2.4 3.4 1.5 2.5 3 . 5 El progranrn sidsiesi raknlil i escriu sobre 1111 fitx Pr el valor de la funció objecti11, gra.clieut i hf'ssiana. solm " el punt. :r coutingut a.l fitxPr sids. var. Una sessió co1npleta (l'ncpwst progn1111a sPri<i: ....... ' .. .. ... . ......... . . . . . . . . . . . . . . . . . . . . . . . . : : ONL, Dept. EIO. UPC. OPTIMITZACIO DE MODELS NO LINEALS SENSE CONSTRICCIONS 'SUDDEN INFANT DEATH SYNDROME (SIDS)' Comprovacio de les rutines d'usuari. >Numero d'identificacio? : JOO(Retu.rnj > Serie que es vol ajustar (morts=1, temp.=2)? 2 (Ret11.rn} > Generacio de dades. > Lectura fitxer sids.var > -- > -=- > :L ., -- -; 7 P<" =- --- · -~ó ·-!..a 8 D1·pt .. ETO / r·¡y· / Op1i111irzacicJ de 111odels 1w lilleals se11..,P coristr iccions > Calcul de f.o. i gradient > Calcul de l'hessiana > ok. El fitxer de resultats es diu sids100.test SORT ! ! ! FORTRAN STOP en aque~t ca s 8'lm demanat <'l ,.,\ kul de f(:c) , 'Vf( :r ) i \7 2 f(.7:) pel prol>k ma número 100 concspouent. a la s i· rie d<" t.empcratur es. El fit.xer sids10 0. test creat és: : : Grup ONL, Dept. EIO . UPC . ...... ... . ... .. . ... ... .. . .. ..... .. .......... ... ... . .. OPTIMITZACIO DE MODELS NO LINEALS SENSE CONSTRICCIONS PRACTICA 'SUDDEN INFANT DEATH SYNDROME (SIDS)' Comprovacio de les rutines d'usuari del problema 100 . : : Punt de prava : ::: : . . ::"::: ::::: x( 1) =0.llOOOOOE+Ol x( 4) = 0.1200000E+Ol x( 7) =0 . 1300000E+Ol x(10)=0.1400000E+Ol x(13)=0.15 000DOE +Ol x( 2)=0.2 100000E+Ol x( 5 )=0 2200000E+01 x( 8) =0 2300 000E+Ol x(ll)=O 2400000E+Ol x(14)=0 250 0000E+01 :::::::::::::: ::::: ........ .. : : Funcio objectiu .. ....... . .... = 0 .8121 322E+05 g( 1)= -. 1141400E+04 g( 4) =-. 1101335E+04 g( 7) =- .1135874E+04 g(10)=-.9297095E +0 3 g(13)=-.8615359E+03 : : Hessiana : : Elements diago nals gradient g( 2)=0 4542151E+03 g( 5)=0 9661762E+03 g( 8)=0.1169051E+04 g(ll) =0.14 15346E+04 g(14)=0 .1806016E+04 x( 3)=0.3100000E+Ol x( 6)=0.3200000E+Ol x( 9)=0.3300000E+Ol x(12 )=0 .3400000E+Ol x(15)=0.3500000E+01 g( 3)=-.6407929E +03 g( 6)=-.8856520E +03 g( 9)=-.9996044E+03 g(12)=-.8431335E+03 g(15)=- .9206823E+03 hd( 1) = 0.1825000E+03 hd( 2) = 0 . 1476365E+04 hd( 3) = 0 .6 679500E+05 hd( 4) = 0. 18 25000E+03 hd( 5) = 0 . 1584402E+04 hd( 6) = 0 .6679500E+05 hd( 7) = 0.182 5000E+03 hd( 8) = 0.1785061E+ 04 hd( 9) = 0.6679500E+OS hd ( 10) = 0. 1825000E+03 hd( 11) = 0.1659293E +04 hd( 12) = 0 .6679 500E+05 hd( 13) = 0.1825000E+03 hd ( 14) = 0.1702929E+04 hd( 1 5) = 0 .66795 00E+05 2 Doc1111u•11l : ació de fa r11ti1Ja Efl. Jf .B F users should pass information to FUNCT via COMMON and not use IW or W at all. In any case FUNCT must not change the first 2 elements of IW or the first max(8,7•N+N•(N-1)/2) elements of w. ••Note••: FUNCT should be tested separately before being used in conjunction with //E04LBF//. It must be declarad as EXTERNAL in the (sub)program from which //E04LBF// is called. Parameters denoted as Input must ••not** be changed by this procedure . 3: HESS - SUBROUTINE, supplied by the user. External Procedure HESS rnust calculate the second derivatives of F at any point x. (As with FUNCT, there is the option of causing //E04L BF// to terrninate inunediately. ) Its specification is: SUBROUTINE HESS(IFLAG, N, XC, FHESL, LH, FHESD, IW, LIW, W, LW) INTEGER IFLAG, N, LH, IW(LIW), LIW, LW //real// XC(N), FHESL(LH), FHESD(N), W(LW) 1: IFLAG - INTEGER. On entry: !FLAG is set to a non-negative number. Input/Outpu t On exit: if HESS resets !FLAG to sorne negative number, // E04LBF// will terminate immediately with IFAIL set to the user's setting of !FLAG. 2: N - INTEGER. Input On entry: the nurnber n of variables. 3: XC(N) - //real// array. Input On entry: the point x at which the second derivatives of F are required. 4: FHESL(L H) - //real// array. Ou tp ut On exit: unless IFLAG is reset, HESS must place the strict lower triangle of the second derivativa rnatrix of F (evaluated at the point x) in FHESL, stored by rows, i.e. set FHESL((i-l)(i-2)/(2+j )) = ((pd••2 F)/((pd x(i)) ( pd x(j)))) at the point XC , far i = 2,3, ... ,n; j = 1,2, ... ,i-1. (The upper triangle is not required because the rnatrix is syrnmetric. ) 5: LH - INTEGER . Input On entry : the len gth of t he array FHESL. 6: F HE S D(N) - //real// arra y. Input/Out put On e ntry: the valu e of (pd F)/(pd x(j)) at the point x, for j = 1,2, .. .,n . These values may be useful in the evaluation of the second derivatives. On exit: unless IFLAG is reset, HESS must place the diagonal el eme nts of the second deri va tive matrix of F ( evaluated at t he point x) in FHE SD , i.e. set FHESD(j) = ((pd••2 F) /( pd x(j)••2)) at the point XC , fo r j = 1 ,2 , ... 1 n. 7: IW(LIW) - INTE G ER a rr ay. 8: LIW - INTEGER. \tlork space Input 15 16 :'._j ,..!" ·! ..:. .:'.!: : t.~ ~:. lG Dept .EIO / l "PC' / Opri111itz;u.·úi tle rnm!els no liw-•als sense co1i strícrio11s 9: W(LW) - //real // array. 10: LW - INTEGER. Workspace Input As in FUNCT , these parameters co r respond to the parameters IW, LIW, W, LW of //E04LBF// . HESS ••must not change•• the first two elements of IW or the first roax(8, 7•N+N•(N-1)/2) elements of W. Again, it is reconunended that the user should pass quantities to HESS via COMMON and not use IW or W at a ll. ••Note ••: HESS should be tested separately before being used in con junction with //E04LBF//. It must be declarad as EXTERNAL in the (sub)prograro from wh ich //E04LBF// i s callad . Pararoete rs denoted as Input rnust **not** be changed by this procedure. 4: MONIT - SUBROUTINE, supplied by the user. Externa l Procedure If IPRINT >= O, the user must supply a subroutine MONIT which is suitable for monit oring the minimization process. MON IT must not change the values of any of its parameters . If IPRINT < O, a ro utine MONIT with the correct pararoeter list should st i ll be supplie d, although it will not be called. Its specitication is: SUBROUTINE MONIT(N, XC, FC , GC , ISTATE, GPJNRM, COND, PO SDEF , 1 NITER, NF, IW , LIW, W, LW) LOGICAL POSDEF INT EGER N, ISTATE(N), NITER, NF, IW (LIW), LIW, LW // real/ / XC(N), FC , GC(N), GPJNRM, COND, W( LW) 1: N - INTEGER. On entry: the number n of variable s. 2: XC(N) - //real // array. On entry: the co-ordinates of the current poin t x. 3: FC - //rea l //. On entry: the v alue of F(x) at the current po int x. Input Input Input 4: GC(N) - //real // arra y. Input On entry : the value of (pd F)/(pd x(j)) at the current point x, for j = 1 ,2 , ... ,n. 5: ISTATE(N) - INTEGER array. Input On entry: informati on about which variables are currently fixed on their bounds and which are free. If ISTATE(j) is ne gativ e, x(j) is currently: - fixed on its u pp er bound if ISTATE(j) = -1 - fi xed on its lo w er bound if ISTATE(j) = -2 - effectively a constan t (i.e. l(j) = u( j) ) if ISTATE(j) = -3. If ISTA TE is positive, it s value g ives the posit ion of x(j) in the sequence of free va riables. 6: GPJNRM - //re al//. Input On entry : the E uc li dean n orrn of the projected gr a dient v ecto r g(z). 7: COND - //real // . Inp ut On entry: the ratio of the largest t o th e smal le st elements 2 Donw1e11t;1ciá dt· la r11ri11a EO-lLBF of the diagonal factor D of the projected Hessian matrix (see specification of HESS below). This quantity is usually a good estimate of the condition number of the projected Hessian matrix. (If no variables are currently free, CONO is set to zero.) 8: POSDEF - LOGICAL. Input On entry : POSDEF is set .TRUE. or .FALSE. according to whether the second derivative matrix for the current subspace, H, is positive-definite or not. 9: NI TER - INTEGER . On entry: the nurnber of iterations (as outlined in of the routine document in the NAG Fortran Library which have been performed by //E04LBF// so far. 10: NF - INTEGER. Input Section 3 Manual) Input On entry: the number of times that FUNCT has been called so far. Thus NF is the number of function and gradient ev aluations made so far. 11: IW(LIW) - INTEGER array. 12: LIW - INTEGER. 13: W(LW) - //real// array. 14: LW - INTEGER. Workspace Input Workspace Input As in FUNCT, and HESS, these parameters correspond to the parameters IW, LIW, W, LW of //E04LBF//. They are included in MONIT's parameter list primarily for when //E04LBF// is called by other library routines. The user should normally print out FC, GPJNRM and CONO so as to be able to compare the quantities mentioned in Section 7 of the routine document in the NAG Fortran Library Manual. It is usually helpful to examine XC, POSDEF and NF as well. MONIT must be declared as EXTERNAL in the (sub)program from which // E04LBF// is called. Parameters denoted as Input must **not•* be changed by this procedure. 5: IPRINT - INTEGER. Input On entry: the frequency with which MONIT is to be called. If IPRINT > O, MONIT is called once every IPRINT iterations and just befare exit from / /E04LBF// . If IPRINT" O, MONIT is just called at the final point. If IPRINT < O, MONIT is not called at all. IPRINT sho uld normally be set to a small positive number. Suggested value: IPRINT " l. 6: MAXCAL - INTEGER. Input On entry: the maximum permitted number of evaluations ot F(x), i.e . the maxirnum perrnitted nurnber of calls of FUNCT . Suggested valu e: MAXCAL" SO•N. Constraint: MAXCAL >; 1. 17 18 :: :...:.._~- _::.i..:=: :: _ f_;:f- - : -=.:: _:- • •:: 1 :'::. ~-- ;:.:::-:i. :: : ::- 17 18 ].!:!: ;;: ...::. - ' ·- l.1"7.:Z..: !. c::·i:: .!.!: _f'_;_::..::. -- -- - ~:: ""-::57 L.-=-=- x :. Depr.EJO / 1 ·p( '/ <J¡>riIJ1itz;1ciú de 1110tlt•ls 110 liueals se11sr~ <'on..;trícc io ns 7: ETA - //real//. Input On entr y: every iteration of //E04LBF// invo lves a linear minimization (i . e. mini mizatio n of F(x+alpha•p) with respect to alpha). ETA specifies how accurately these linear rninimizations are to be performed. The mínimum with respect to alpha will be locat ed more accurately for srnall values of ETA ( say 0.01) than for large values (say 0.9). Although accurate linear min imizations wil l generally reduce the nurnber of iterations of //E04LBF//, this usually results in an increase in the number of func tion and gradient evaluations required for each iterat io n. On balance, it is usually more efficient to p erform a low accuracy line ar minimization . Suggested value: **ETA ; 0.9 is usually a good choice** although a smaller value may be warranted if the rnatrix of second derivatives is expensive to compute comparad with the function and first derivatives. ••If N = !, ETA should be set to O.O•• (also when the problern is effectively 1-dimens ional even though n > 1; i.e . it for all except one of the variables the lower and upper bounds are equal). Constrai n t: O.O<= ETA< ! . O. 8: XT OL - //real//. Input On entry: the accuracy in x to whi ch the solution is required. If x(true) is the true value of x at the mínimum, then x(sol), the estimated position prior to a normal exit, is such that (norrn of x(sol)-x(true) < XTOL•(!.O+norrn of x(true))), where (norrn of y= SQRT(SUM frorn j=l to n((y(j))••2))). For example, if the eleme n ts of x(sol) are not much larger than !.O in rnodulus, and i f XTOL is set to !0••(-5) then x(sol) is usually accurate to about 5 decimal places. (For further details see Section 7 of t he routine document in the NAG Fortran Lib raryManual.) If the problem is scaled roughly as described in Section 8 of the routine docurnent in the NAG F ort ran Library Manual and epsilon is the //rnachine precision//, then SQRT(epsilon) is probably the srnallest r easo nable choice for XTOL. (This is because, normally , to machine accuracy, F(x+(SQRT(epsilon)) e(J)) = F(x) wh ere e(j) is any colurnn of th e iden tity m atrix.) If the user sets XTOL to O.O ( or any positive value le ss than epsil on), //E04 LBF // will use 10.0•(SQRT(epsilon)) inste ad of XTOL. Sugg ested valu e: XTOL = O.O. Constraint: XTOL >= O.O. 9 : STEPMX - // real/ /. Input On entry: an estímate of the Euclidean distance between the solution and the startin g poi nt supplied by the us e r. (For maximurn efficiency a sl igh t overest imat e is pre fera ble.) / / E0 4LBF// will en sure that, for each i teration, SQRT(SU M fro rn j =l to n(( (x(J)(k)) - (x(j)(k- 1 ))) •• 2) ) <= STEPMX 2 Doci1111f'11tarití ch" Ja r11tí11a EUJ!. HF where k is the iteration num b er. Thus, if the problem has more than one solution, //E04LBF// is rnost likely to find the one nearest to the starting point. On difficult problems, a realistic choice can prevent the sequence of x(k) entering a region where the problem is ill-behaved and can also help to avoid possible overflow in the evaluation of F(x). However an underestimate of STEPMX can lead to inefficiency. Suggested value: STEPMX = 100000.0. Constraint: STEPMX >= XTOL. 10: IBOUND - INTEGER. Input On entry: specifies whether the problem is unconstrained or bounded. If there are bounds on the variables, IBOUND can be used to indicate whether the facility for dealing with bounds of special forros is to be used. IBOUND should be set to one of the following values: IBOUND = O if the variables are bounded and the user will be supplying all the l(j) and u(j) individually. IBOUND = 1 if the problem is unconstrained. IBOUND = 2 if the variables are bounded, but all the bounds are of the forrn O<= x(j). IBOUND = 3 if all the variables are bounded, and 1(1) = 1(2) = = l(n) and u(l) = u(2) u(n). IBOUND = 4 if the problern is unconstrained. (The IBOUND 4 option is provided purely for consistency with other routines. In //E04LBF// it produces the sarne effect as IBOUND = 1.) Constraint: O<= IBOUND <= 4. 11: BL(N) - //real// array. Input/Output On entry: the fixed lower bounds l(j). If IBOUND is set to O, the user rnust set BL(j) to l(j), for j = 1,2, ... ,n. (If a lower bound is not specified for any x(j), the corresponding BL(j) should be set to a large negativa nurnber, e.g. -10••6.) If IBOUND is set to 3, the user rnust set BL(l) to 1(1); //E04LBF// will then set the rernaining elernents of BL equal to BL(l). If IBOUND is set to 1, 2 or 4, BL will be initialised by //E04LBF//. On exit: the lower bounds actually used by //E04LBF//, e.g. If IBOUND = 2, BL(l) = BL(2) = BL(N) = O.O. 12: BU(N) - //real// array. Input/Output On entry: the fixed upper bounds u(¡). If IBOUND is set to O, the user rnust set BU(j) to u(j), for j = 1,2, ... ,n. (If an upper bound is not specified for any variable, the corresponding BU(j) should be set to a large positive number, e.g. 10**6.) 20 _;:. ., ~ -- -' : ~ -' ., ..::. .. ! - -- -• :!J _•J -- ~ ~ - ~: . = - ~:! -=. .. - : ~ - 1 - 1 - 1 - l [Q 20 r ~ Dt•¡Jt.F.10 / r·p c / Opti111i1z;u· i<J dt:> lJHHlt~Js llO liJJeals St'/JSe co11striccions If IBOUND is set to 3, the user must set BU(!) to u(l); //E04LBF// will then set the remaining elements of BU equal to BU(!). If IBOUND is set to 1, 2, or 4, BU will then be initialised by //E04LBF//. Dn exit: the upper bounds actually used by //E04LBF//, e.g. If IBOUND = 2, BU(!) = BU(2) = = BU(N) = 10••6. 13: X(N) - //real// array. Input/Output On entry: X(j) must be set to a guess at the jth component of the position of minimum, for j = 1,2 ... . ,n. Dn exit: the final point x(k). Thus, if IFAIL =O on exit, X(j) is the jth component of the estimated position of the rninirnum. 14: HESL(LH) - //real// array. Output See description of HESD below. 15: LH - INTEGER. Dn entry: the actual length of HESL as declarad in the (sub)program from which //E04LBF// is called. Constraint: LH >= max(N•(N-1)/2,1). Input 16: HESD(N) - //real// array . Output On exit: during the determination of a direction p(z) (see Section 3 of the routine document in the NAG Fortran Library Manual), H +E is decomposed into the product L•D•(transpose(L)), where Lis a unit lower triangular rnatrix and D is a diagonal rnatrix. (The matrices H, E, L and O are all of dirnension n(z). where n(z) is the number of variables free from their bounds. H consists of those rows and columns of the full second derivative matrix which relate to free variables. E is chosen so that H +E is positive-definite.) HESL and HESD are used to store the factors L and D. The elements of the strict lower triangle of L are stored row by row in the first n(z)(n(z ) -1 )/2 positions of HESL. The diagonal elements of D are stored in the first n(z) positions of HESD. In the last factorization before a normal exit , the matrix E will be zero, so that HESL and HESD will contain on exit, the factors of the final second derivative matrix H. The elements of HESD are useful for deciding whether to accept the result produced by //E04LBF// (see Section 7 of the routine document in the NAG Fortran Library Manual). 17: ISTATE(N) - INTEGER array . Output Dn exit: information about which variables are currently on their bounds and which are free. If ISTATE(j) is: - equal to -1, x(j) is fixed on its upper bound - equal to -2, x(j) is fixed on its lower bound - equal to -3, x(j) is effectively a constant (i.e. l(j) = u(j)) - positive, ISTATE(j) g ive s the position of x(j) in the sequence of free variables. 2 Docr1111e1daci1) dt• la r11ti11;1 ELJ .Jl.13F 18: F - //real//. Output On exit: the function value at the final point given in X. 19: G(N) - //real// array. Output On exit: the first derivative vector corresponding to the final point given in X. The components of G corresponding to free variables should normally be clase to zero. 20: IW(LIW) - INTEGER array. Workspace 21: LIW - INTEGER. Input On entry: the length of IW as declared in the (sub)program from which //E04LBF// is callad. Con s traint: LIW >= 2. 22: W(LW) - //real// array. 23: LW - INTEGER. On entry: the length of W as declarad in which //E04LBF// is called. Constraint: LW >= max(7•N+N• (N -l)/2,8). 24: IFAIL - INTEGER. Workspace Input the (sub)program from Input/Output On entry: IFAIL must be set to O, - 1 or 1. U sers who are unfamiliar with this parameter should refer to Chapter P01 in the NAG Fortran Library Manual or this HELP system for details. On exit: IFAIL ; O unless the routine detects an error or gives a warning (see Section D). **For this routine**• because the values of output parameters may be useful even if IFAIL<>O on exit, users are recom.mended to set IFAIL to -1 before entry. ••It is then essential to test the value of IFAIL on exit**· To suppress the output of an error rnessage when soft failure occurs , set IFAIL to 1. 21 22 = = = = = '.?l ~ !': 22 Dept..EIO /(·pe/ O¡>ti11iitzacitJ de rno dels 1w li11eals se1Jse constriccions 2.1.4 Enor lndicators and Warnings Errors or w ar nings specified by the routine: IFAIL < O A negative value of IFAIL indicates an exit from //E04LBF// because the user has set !FLAG negativa in FUNCT or HESS. The value of IFAIL will be the sama as the user's setting of IFLAG . IFAIL = 1 On entry, N < 1, or MAXCAL < 1, or ETA < O.O, or or or or or or or b ETA >= 1. O, XTOL < O O, STEPMX < XTOL, IBOUND < O, IBOUND > 4, BL(j) > BU (j) fo r sorne j if IBOUND BL(l) > BU (!) if IBOUND = 3 , or LH < rnax(!,N•(N-1)/2), or LIW < 2, or LW < rn ax(8,7•N+N•(N- 1) /2). o, (Note that if the u ser has set XTOL to O.O, //E0 4LBF// u ses the default value and continues without failing.) When this exit occurs no values will have been assigned to F or to the elernents of HESL, HESD or G. IFAIL = 2 There have been MAXCAL function evaluations. If steady reducti o ns in F(x) were monitored up to the poi nt where this exit occurred, then the exit probably occurred simply because MAXCAL was set too small, so the calculations should be restarted from th e final point held in x. This exit rnay also indicate that F(x) has no rninimum. IFAIL = 3 The conditions for a minimum have not all been rnet 1 but a lover point could not be found. P rovided that, on exit, the first derivativas of F(x) with respect to the free variables are suff icien t ly small, and that the estimated condition number of the second derivativa rnatrix is not too large , this error exit rnay simply mean that, although it has not been possible to satisfy the s pecified requirements, the algorithrn has in fact found th e rninimum as far as the accuracy of the machine perrnits. Such a s ituati on can arise, for instance, if XTOL has been set so small that rounding errors in the evaluation of F(x ) or its derivatives make it irnpossible to satisfy the convergence condit ions. If the estirnated condition number of the second derivative matrix at the final point is l arge , it could be that the final point is a mínimum, but that the smal lest eigenvalue of th e Hessi an m at rix is so clase to zero that it is not possible to recognise the point as a mínimum. 2 Docrw1<~utaciú dt.• la r11ti1Ja EIJ.ILBF IFAIL = 4 !lot used. (This is done to make the signíficance of IFAIL = 5 similar for //E04JBF//, //E04KBF//, //E04KDF// and //E04LBF//.) IFAIL = 5 All the Lagrange-multiplier estimates which are not indisputably positive lie relatively close to zero, but it is irnpossible either to continua minimizing on the current subspace or to find a feasible lower point by releasing and perturbing any of the fixed variables. The user should i nvestigate as for IFAIL = 3. The values IFAIL = 2, 3 and 5 may also be caused by mistakes in FUNCT or HESS, by the formulat1on of the problem or by an awkward function. If there are no such mistakes, it is worth restarting the calcu lati ons from a diff erent starting point (n ot the point at which the failure occurred) in arder to avoid the regían which caused the failure. END OF E04LBF FORTRA N SUMMARY - MARK 15 MARCH 1991 2.2 Example Program Source Code 2.2 .1 Main prograrn E04LBF Example Program Text. Mark 14 Revised. NAG Copyright 1989. Parameters INTEGER PARAMETER INTEGER PARAMETER . . Local Scalars N, LH, LIW, LW (N=4,LH=N•(N-1)/2,LI\1=2,L\1=7•N+N•(N-1)/2) NOUT (NOUT=6) DOUBLE PRECI SI DN ETA , F, STEPMX, XTOL INTEGER IBOUND, IFAIL, IPRINT, J, MAXCAL Local Arrays .. OOUBLE PRECISION BL(N), BU(N), G(N), HESD(N), HES L(LH), W(LW), + X(N) INTEGER ISTATE(N), IW(LIW) External Subroutines EXTERNAL E04HCF, E04HDF, E04LBF , FUNCT, HESS, MONIT . . Executable Statements .. WRITE (NOUT,•) 'E04LBF Example Program Re sults ' '.?4 23 _J 24 DepLEIO / (. ' }.J (' / O¡ >tirnitzaciú de models 110 lineals se11se co11s tric c ions Set up X(1) X(2) X(3) X(4) an arbitrary point at which to check the derivativas 1.4600 -0 .8200 0.5700 1.2100 Check the lst derivatives <3> IFAIL = O CALL E04HCF(N,FUNCT,X,F,G,IW,LIW,W,LW,IFAIL) Check the 2nd derivatives <4> IFAIL = O CALL E04HOF(N,FUNCT,HESS,X,G,HESL,LH,HESO,IW,LIW,W,LW,IFAIL) Continua setting parameters f or E04LBF * Set IPRINT to 1 to obtain output frorn MONIT at each iteration * IPRINT = -1 MAXCAL = SO•N ETA = 0.900 Set XTOL to zero so that E04LBF will use the default tolerance XTOL = 0.000 We estimate that the minimum will be within 4 units of the starting point STEPMX = 4.000 IBOUNO = O BL ( 1) = 1. 000 <r» BU(1) = 3.000 BL(2) = -2.000 BU(2) = 0.000 X(3) is not bounded, so we set BL(3) to a larga negativa nurnber and BU(3) to a large positiva number BL(3) = -1. 006 BU(3) = 1. 006 BL ( 4) = 1. 000 BU(4) = 3.000 Set up X(1) X(2) X(3) X(4) starting point 3 .000 -1.000 0.000 1.000 IFAIL = 1 <J> En el nost.re ca.-;, aq1ws1 a co111prm·;1riú no la fart"' lll <4> ídem <f;> A part.ir d' aq11í , fi11 s a .. Set up starting point " s' i11lrialit .z t>n fit.P- s. Com <JtlP ni (Pt) ni ( Pm) en t e tH -' ll, aquesta par t la podt-'111 eliiuinar