donlp2.c

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

        acond = epsmac/tolmac;
    }
    den1 = pow(o8vecn(1,n,ltdx),2);
    den2 = o8sc1(1,n,dg,difx);
    if ( den1 <= rho1*anorm*pow(ndx,2) || acond >= one/rho1 ) {
    
        /* take a restart step */
        
        o8inim();
        
        return;
    }
    if ( nres == 0 ) {
    
        /* in the unconstrained case we take the Powell update */
        
        th = one;
        if ( den2 < p2*den1 ) {
            th = p8*den1/(den1-den2);
            for (i = 1 ; i <= n ; i++) {
                dg[i] = th*dg[i]+(one-th)*adx[i];
            }
            den2 = o8sc1(1,n,dg,difx);
        }
        term = one/sqrt(den2);
        for (i = 1 ; i <= n ; i++) {
            dg[i]   = dg[i]*term;
            updz[i] = dg[i];
        }
        term = one/sqrt(den1);
        for (i = 1 ; i <= n ; i++) {
            updx[i] = adx[i]*term;
        }
        accinf[itstep][28] = den2/den1;
        accinf[itstep][29] = th;
        accinf[itstep][27] = two;
        if ( th != one ) accinf[itstep][27] = three;
    } else {
        ngtdx = o8vecn(1,alist[0],gtdx);
        term  = one/sqrt(den1);
        for (i = 1 ; i <= n ; i++) {
            updx[i] = adx[i]*term;
        }
        if ( den2 >= rho1*o8sc1(1,n,dg,dg)
        && o8vecn(1,n,dg) >= sqrt(epsmac)*ndx ) {
            xsik = zero;
            for (i = 1 ; i <= n ; i++) {
                updz[i] = dg[i];
            }
            den21 = den2;
        } else {
        
            /* try Pantoja-Mayne modification */
            
            den3 = tk*pow(ndx,2)+pow(ngtdx,2);
            if ( den2 >= rho1*o8sc1(1,n,dg,dg) ) {
                xsik = one;
            } else {
                xsik = one+(tk*pow(ndx,2)+fabs(den2))/den3;
            }
            for (i = 1 ; i <= n ; i++) {
                term = zero;
                for (j = 1 ; j <= alist[0] ; j++) {
                    term1 = gres[i][alist[j]]*gtdx[j];
                    term1 = term1/gresn[alist[j]];
                    term  = term+term1;
                }
                updz[i] = dg[i]+xsik*(tk*difx[i]+term);
            }
            den21 = o8sc1(1,n,updz,difx);
        }
        term = one/sqrt(den21);
        for (i = 1 ; i <= n ; i++) {
            updz[i] = updz[i]*term;
        }
        th = one;
        if ( den2 < p2*den1 ) {
            th = p8*den1/(den1-den2);
            for (i = 1 ; i <= n ; i++) {
                dg[i] = th*dg[i]+(one-th)*adx[i];
            }
            den2 = o8sc1(1,n,dg,difx);
        }
        term = one/sqrt(den2);
        for (i = 1 ; i <= n ; i++) {
            dg[i] = dg[i]*term;
        }
        if ( o8vecn(1,n,dg) <= tm3*o8vecn(1,n,updz) ) {
        
            /* the Powell update produces a smaller growth */
            
            for (i = 1 ; i <= n ; i++) {
                updz[i] = dg[i];
            }
            accinf[itstep][28] = den2/den1;
            accinf[itstep][29] = th;
            accinf[itstep][27] = two;
            if ( th != one ) accinf[itstep][27] = three;
        } else {
        
            /* no update if strongly irregular */
            
            accinf[itstep][27] = one;
            accinf[itstep][28] = tk;
            accinf[itstep][29] = xsik;
        }
    }
    o8upd(a,updz,updx,n,&fail);
    
    /* check illconditioning after updating */
    
    term  = fabs(a[1][1]);
    term1 = term;
    i     = 1;
    
    /* in order to overcome a curious error in hp's f77 compiler */
    /* this kind of loop                                         */
    
    while ( i < n ) {
        i     = i+1;
        term  = max(term, fabs(a[i][i]));
        term1 = min(term1,fabs(a[i][i]));
    }
    if ( fail || pow(term1,2) <= rho1*pow(term,2) ) {
    
        /* reset */
        
        o8inim();
    }
    return;
}

/* **************************************************************************** */
/*                      write short information on standard out                 */
/* **************************************************************************** */
void o8shms(void) {

    static DOUBLE   umin;

    if ( te0 && ! silent ) {
        umin = accinf[itstep][11];
        printf(
        "%5i fx= %14.6e upsi= %8.1e b2n= %8.1e umi= %8.1e nr%4i si%2i\n",
        itstep,fx,upsi,b2n,umin,nr,(int)accinf[itstep][10]);
        
        fprintf(prou,
        "%5i fx= %14.6e upsi= %8.1e b2n= %8.1e umi= %8.1e nr%4i si%2i\n",
        itstep,fx,upsi,b2n,umin,nr,(int)accinf[itstep][10]);
    }
    return;
}

/* **************************************************************************** */
/*                      write messages on "special events"                      */
/* **************************************************************************** */
void o8msg(INTEGER num) {

    static INTEGER  i;

    if ( num <= 0 || num  > 25 ) return;
    
    switch (num) {
    
    case 1:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"rankdeficiency of grad's of active constr.!\n");
        fprintf(meu,"on the basis of delmin!\n");
        
        return;
        
    case 2:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"rescaling of objective function scf= %.15e\n",scf);
        
        return;
        
    case 3:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"rankdeficiency of grad's of active constr.!\n");
        fprintf(meu," del= %.15e\n", del);
        
        return;
        
    case 4:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"rankdeficiency of grad's of active constr.!\n");
        fprintf(meu," del= %.15e\n", del);
        
        return;
        
    case 5:
        fprintf(meu,"qpterm<0: no dir. of. desc., tauqp max\n");
        
        return;
        
    case 6:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"second order correction suppressed! \n");
        
        return;
        
    case 7:

        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"retry next step with a=id \n");
        
        return;
        
    case 8:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"some constraint has gradient equal to zero \n");
        fprintf(meu,"resulting d may be no direction of descent\n");
        
        return;
        
    case 9:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"try regularized SQP with increased weights\n");
        fprintf(meu,"since dnorm small or infeasibility large\n");
        
        return;
        
    case 10:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"QPsolver did not complete, try d anyway, may fail\n");
        
        return;
        
    case 11:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"direct. deriv. positive! may be due to inaccurate\n");
        fprintf(meu,"gradients or extreme illconditioning\n");
        fprintf(meu,"dirder=  %.15e\n", dirder);
        
        return;
        
    case 12:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"call of matdru suppressed, mat too large\n");
        
        return;
        
    case 13:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"startvalue corrected in order to fit bounds\n");
        
        return;
        
    case 14:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"try full SQP due to slow progress in x \n");
        
        return;
        
    case 15:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"try full SQP due to small stepsizes while\n");
        fprintf(meu,"infeasibility large\n");
    case 16:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"on exit from o8qpdu dir. deriv. positive!\n");
        fprintf(meu,"try increased tauqp\n");
        
        return;
        
    case 17:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"try regularized SQP with increased weights\n");
        fprintf(meu,"no decrease of weights possible\n");
        fprintf(meu,"and incompatibility large\n");
        
        return;
        
    case 18:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"try regularized SQP with increased weights\n");
        fprintf(meu,"since no direction of descent has been obtained\n");
        
        return;
        
    case 19:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"degeneracy in dual QP\n");
        fprintf(meu,"restr. %i to be added\n",iptr);
        fprintf(meu,"new rii= %.15e\n",riitr);
        fprintf(meu,"length of current working set=%i\n",iqtr);
        fprintf(meu,"working set\n");
        for (i = 1 ; i <= iqtr ; i++) {
            fprintf(meu,"%4i ",aitr[i]);
            if ( i % 15 == 0 || i == iqtr ) fprintf(meu,"\n");
        }
        fprintf(meu,"primal feasibility violation is= %.15e\n",sstr);
        
        return;
        
    case 20:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"dual QP seemingly infeasible \n");
        fprintf(meu,"theoretically impossible\n");
        
        return;
        
    case 21:
        fprintf(meu,"step=%i\n", itstep);
        fprintf(meu,"no update since dg=0\n");
        
        return;
        
    case 22:
        fprintf(meu,"step%i\n", itstep);
        fprintf(meu,"function evaluation returns err=true\n");
        
        return;
        
    case 23:
        fprintf(meu,"step%i\n", itstep);
        fprintf(meu,"contraint evaluation returns err=true\n");
        
        return;
        
    case 24:
        fprintf(meu,"step%i\n", itstep);
        fprintf(meu,"current point cannot be changed in current\n");
        fprintf(meu,"direction due to error-message from function\n");
        fprintf(meu,"evaluation\n");
        
        return;
        
    case 25:
        fprintf(meu,"step%i\n", itstep);
        fprintf(meu,"reduce stepsize due to error-indicator set\n");
        fprintf(meu,"by users function evaluation\n");
        
        return;
    }
}

/* **************************************************************************** */
/* equality constrained recursive quadratic programming with                    */
/* multiple inactivation and superlinearly convergent projected                 */
/* BFGS-update (version 12/93 Spellucci )                                       */
/* **************************************************************************** */
void o8opti(void) {

    void   o8info (INTEGER icase);
    void    o8sce  (void);
    void    o8bfgs (void);
    void    o8shms (void);
    void    o8msg  (INTEGER num);
    void    o8inim (void);
    void    o8dird (void);
    void    o8cutd (void);
    void    o8smax (void);
    void    o8unim (DOUBLE sig1th);
    void    o8egph (DOUBLE gphi[]);
    void    o8dec  (INTEGER nlow,INTEGER nrl);
    void    o8ht   (INTEGER id,INTEGER incr,INTEGER is1,INTEGER is2,INTEGER m,
                    DOUBLE a[][NRESM+1],DOUBLE beta[],DOUBLE b[],DOUBLE c[]);   
    void    o8sol  (INTEGER nlow,INTEGER nup,DOUBLE b[],DOUBLE x[]);
    void    o8solt (INTEGER nlow,INTEGER nup,DOUBLE b[],DOUBLE x[]);
    void    o8rght (DOUBLE a[][NX+1],DOUBLE b[],DOUBLE y[],DOUBLE *yl,INTEGER n);
    void    o8left (DOUBLE a[][NX+1],DOUBLE b[],DOUBLE y[],DOUBLE *yl,INTEGER n);
    DOUBLE  o8vecn (INTEGER nl,INTEGER nm