donlp2.c

sqp程序包。用sqp算法实现非线性约束的优化求解

/* Conditions of use:                                                           */
/* 1. donlp2 is under the exclusive copyright of P. Spellucci                   */
/*    (e-mail:[e-mail])                             */
/*    "donlp2" is a reserved name                                               */
/* 2. donlp2 and its constituent parts come with no warranty, whether ex-       */
/*    pressed or implied, that it is free of errors or suitable for any         */
/*    specific purpose.                                                         */
/*    It must not be used to solve any problem, whose incorrect solution        */
/*    could result in injury to a person , institution or property.             */
/*    It is at the users own risk to use donlp2 or parts of it and the          */
/*    author disclaims all liability for such use.                              */
/* 3. donlp2 is distributed "as is". In particular, no maintenance, support     */
/*    or trouble-shooting or subsequent upgrade is implied.                     */
/* 4. The use of donlp2 must be acknowledged, in any publication which contains */
/*    results obtained with it or parts of it. Citation of the authors name     */
/*    and netlib-source is suitable.                                            */
/* 5. The free use of donlp2 and parts of it is restricted for research purposes*/
/*    commercial uses require permission and licensing from P. Spellucci.       */

/* d o n l p  2                                                                 */

/*    version 29/05/98 (*)                                                      */
/*    tauqp dependent on scalres only                                           */
/*    weights computed in a modified version in the singular case               */
/*    some comparisons relative to zero changed                                 */
/*    error-return for function evaluation as added feature                     */
/*    termination of QP-solver changed (d not set to zero)                      */
/*    new version of BFGS: if nr  =  0 take Powell's update                     */
/*    no suppression of update beforehand (with exception dg = 0)               */
/*    plus some minor corrections                                               */

/*    for consistency reasons variable names cf and cgf changed into            */
/*    icf and icgf                                                              */

/*    added feature numerical differentiation of order 1,2,6                    */
/*    requires new parameters epsfcn, taubnd, difftype                          */

/*    added feature of external blockwise evaluation of functions               */
/*    rather than individual evaluation                                         */
/*    requires new parameter bloc                                               */

/*    new feature of user-supplied scaling of x                                 */

/* **************************************************************************** */
/* (*) direct translation of -, from f77 to ANSI C language                     */
/*     by S. Schoeffert, ASB, Bourges, France                                   */
/* version 14.5.2003 : correction for phase=-1: gfn, gradf undefined removed    */
#include "o8gene.h"

void donlp2(void) {

    void    o8st   (void);
    void    o8fin  (void);
    void    o8opti (void);
    DOUBLE  o8vecn (INTEGER nl,INTEGER nm,DOUBLE x[]);
    void    esgradf(DOUBLE x[],DOUBLE gradf[]);
    void    esgradh(INTEGER i,DOUBLE x[],DOUBLE gradhi[]);
    void    esgradg(INTEGER i,DOUBLE x[],DOUBLE gradgi[]);
    
    void    setup0 (void);
    void    setup  (void);
    
    static DOUBLE   yy[NX+1];
    static INTEGER  i,j;
    
    /* default settings of new parameters */
    
    bloc     = FALSE;
    analyt   = TRUE;
    valid    = FALSE;
    epsfcn   = 1.e-16;
    difftype = 3;
    taubnd   = 1.;
    for (i = 1 ; i <= NX ; i++) {
        xsc[i] = one;
        xtr[i] = zero;
    }
    epsdif = tm8;
    for (i = 0 ; i <= NRESM ; i++) {
        gconst[i] = FALSE;
        val[i]    = FALSE;
        if ( i > 0 ) gresn[i] = one;
    }
    for (i = 1 ; i <= NRESM ; i++) {
        cfuerr[i] = FALSE;
        
        /* the error indicator of the function evaluation */
    }
    ffuerr = FALSE;

    /* some standard initialization which may be overwritten by */
    /* setup0 or setup                                          */

    inx    = FALSE;
    silent = FALSE;
    
    /* the interactive input feature is no longer supported here. for the sake */
    /* of easy revision the variable is supplied here however                  */
    /* if intakt is set TRUE, output to protocol file is copied to stdout in   */
    /* addition                                                                */
    
    intakt = FALSE;
    te0    = FALSE;
    te1    = FALSE;
    te2    = FALSE;
    te3    = FALSE;
    cold   = TRUE;
    
    /* setup0 must initialize analyt, epsdif, del0, tau0 , n , nh , ng , gunit */
    /* gconst , epsfcn , taubnd , analyt , bloc , difftype                     */
    /* and the initial value for x                                             */
    /* may also change settings of all variables initialized above             */

    setup0();

    for (i = 1 ; i <= n ; i++) {
        xst[i] = x[i];
        if ( xsc[i] == zero ) {
            fprintf(stderr,"scaling variable %i is zero\n",i);
            exit(1);
        } 
        x[i] = x[i]/xsc[i];
        yy[i] = xsc[i] ;
    }
    nreset = n;

    o8st();
    
    /* setup may change standard settings of parameters  */
    /* and add some computations in the user environment */

    for (i = 1 ; i <= n ; i++) {
    
        /* ug and og have been evaluted for the original variables */
        /* here we use them internally for the scaled ones         */
        
        if ( llow[i] ) ug[i] = ug[i]/xsc[i];
        if ( lup[i]  ) og[i] = og[i]/xsc[i];
    }
    setup();

    if ( taubnd <= 0 ) {
            fprintf(stderr,"taubnd must not be le zero");
            exit(1);
    }   
    for ( i = 1 ; i <= n ; i++ ) {
       if ( yy[i] != xsc[i] ) {
          fprintf(stderr,"setup has changed xsc, not allowed");
          exit(1);
       }
    }
    
    /* preevaluation of gradients of linear functions */
    /* done only once                                 */

    for (i = 0 ; i <= nres ; i++) {
        if ( gunit[1][i] != 1 && gconst[i] ) {
        
            /* evaluate gradient once */
            
            if ( i == 0 ) {
                val[0] = TRUE;
                
                esgradf(x,gradf);
                
            } else {
                val[i] = TRUE;
                if ( i <= nh ) {
                
                    esgradh(i,x,yy);
                    
                } else {
                
                    esgradg(i-nh,x,yy);
                }
                for (j = 1 ; j <= n ; j++) {
                    gres[j][i] = yy[j];
                }
                gresn[i] = max(one,o8vecn(1,n,yy));
            }
        }
    }
    runtim = clock();
    
    /* call the optimizer */
    
    o8opti();
    
    runtim = (clock()-runtim)/CLOCKS_PER_SEC;

    /* do final solution check and output */
    
    o8fin();
    
    return;
}

/* **************************************************************************** */
/*        initialization program , standard parameter settings done here        */
/* **************************************************************************** */
void o8st(void) {

    void o8msg    (INTEGER num);
    void esh      (INTEGER i,DOUBLE x[],DOUBLE *hxi);
    void esg      (INTEGER i,DOUBLE x[],DOUBLE *gxi);
    void user_eval(DOUBLE xvar[],INTEGER mode);

    static INTEGER  i,j,k,iz;
    static DOUBLE   tol1 ,xnull[NX+1],bd0,infiny,gxi,hxi,term;
    static char     fil[13],xxx[9] = "xxxxxxxx",name1[9];
    static time_t   tim;
    
    epsmac = pow(two,-20);
    
    L100:

    epsmac = epsmac/two;
    term   = one+epsmac;
    
    if ( term != one ) goto L100;
    
    epsmac = epsmac+epsmac;
    tolmac = epsmac;
    
    L200:

    tol1   = tolmac;
    tolmac = tolmac/twop4;
    
    if ( tolmac != zero ) goto L200;
    
    tolmac = tol1;
    
    /* epsmac machine precision, tolmac smallest machine number */
    /* larger than zero (approximately , base 16 for exponent   */
    /* therefore division by 16 assumed)                        */

    /*                        ***** warning *****                        */
    /* on some machines the computation of tolmac may result in an error */
    /* because underflow is not accepted as zero as is assumed here      */

    if ( n > NX ) {
        fprintf(stderr,"donlp2: n too large/recompile\n");
        exit(1);
    }
    if ( nh+ng > NRESM ) {
        fprintf(stderr,"donlp2:nh or ng too large/recompile\n");
        exit(1);
    }
    if ( tau0 == zero ) tau0 = one;
    if ( del0 == zero ) del0 = tau0*p5;
    
    /* if ( del0 > tau0 ) del0 = tau0*p5 */
    
    if ( nreset > n  ) nreset = n;
    if ( nreset <= 4 ) nreset = 4;
    
    /* standard initialization */
    
    lastch = 0;
    lastdw = 0;
    lastup = 0;
    level  = one;
    tau    = tm1;
    iterma = MAXIT;
    epsx   = tm5;
    sigsm  = sqrt(epsmac);
    smalld = tm1;
    
    /* formerly tm2. smalld has much influence on the Maratos-effect */
    
    smallw = exp(two*log(epsmac)/three);
    rho    = tm6;
    rho1   = tm10;
    del01  = del0/tp1;
    delmin = min(del01,max(tm6*del0,smallw));
    if ( ! analyt ) delmin = min(del01,max(epsdif,delmin));
    c1d    = tm2;
    scfmax = tp4;
    taufac = tp1;
    taumax = pow(scfmax,2);
    updmy0 = tm1;
    
    /* take del0 and tau0 from block data or setup0 in function definition */
    /* may be modified by subsequent call of setup                         */
    
    j = 0;
    while ( name[j] == ' ' ) {
        j = j+1;
    }
    if ( name[j] == '\0' ) {
        strcpy(name1,xxx);
    } else {
        k = j+1;
        while ( name[k] != ' ' && name[k] != '\0' && k-j < 8 ) {
            k = k+1;
        }
        strncpy(name1,&name[j],k-j);
        name1[k-j] = '\0';
        
        for (i = 0 ; i <= k-j-1 ; i++) {
            iz = name1[i];
            if ( iz < 48 || ( iz > 57 && iz < 65 )
            || ( iz > 90 && iz < 97 ) || iz > 122 ) name1[i] = 'x';
        }
        if ( k - j < 8 ) strncat(name1,xxx,8-k+j);
    }
    if ( ! silent ) {
        strcpy(fil,name1);
        meu  = fopen(strcat(fil,".mes"),"w");
        strcpy(fil,name1);
        prou = fopen(strcat(fil,".pro"),"w");
    }
    
    infiny = epsmac/tolmac;
    fx     = zero;
    b2n    = zero;
    b2n0   = zero;
    nres   = ng+nh;
    if ( cold ) {
        for (i = 1 ; i <= NX ; i++) {
            for (j = 1 ; j <= NX ; j++) {
                a[j][i] = zero;
            }
            a[i][i]  = one;
            diag0[i] = one;
        }
    }
    for (i = 1 ; i <= NRESM ; i++) {
        diag[i] = zero;
        for (j = 1 ; j <= NX ; j++) {
            qr[j][i]   = zero;
            gres[j][i] = zero;
        }
    }
    for (i = 1 ; i <= NX ; i++) {
        xnull[i] = zero;
        ug[i]    = -infiny;
        og[i]    = infiny;
        llow[i]  = FALSE;
        lup[i]   = FALSE;
    }
    for (i = 1 ; i <= nh ; i++) {
        delfac[i] = one;
    }
    
    if ( bloc ) user_eval(xnull,0);
    
    for (i = nh+1 ; i <= nres ; i++) {
        delfac[i] = one;
        
        /* scan for real lower or upper bounds */
        
        if ( gunit[1][i] == 1 ) {
        
            esg(i-nh,xnull,&gxi);
            
            if ( gunit[3][i] > 0 ) {
                llow[gunit[2][i]] = TRUE;
                ug[gunit[2][i]]   = -gxi/gunit[3][i];
            } else {
                og[gunit[2][i]]   = -gxi/gunit[3][i];
                lup[gunit[2][i]]  = TRUE;
            }