donlp2.c

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

        if ( intakt ) {
            printf(  "\n\n end unimin\n");
            printf(    "\n sig= %11.4e  num. f-evaluations%2i\n", sig,cfincr);
            printf(      " list of inactive hit constraints\n");
            for (i = 1 ; i <= violis[0] ; i++) {
                printf(  "%4i  ", violis[i]);
                if ( i % 13 == 0 || i == violis[0] ) printf(  "\n");
            }
            if ( violis[0] == 0 ) printf("none\n");
        }
        return;
        
    case 11:
        if ( intakt ) printf("additional increase of eta due to large clow\n");
        fprintf(prou,        "additional increase of eta due to large clow\n");
        
        return;
        
    case 12:
        fprintf(prou,
        "\n\n  current scaling,  scf =  %11.4e clow = %12i eta =  %11.4e\n", 
        scf,clow,eta);
        
        if(nres != 0) {
            fprintf(prou,"\n\n  scalres=\n");
            for (i = 1 ; i <= nres ; i++) {
                fprintf(prou,"  %11.4e", w[i]);     
                if ( i % 6 == 0 || i == nres ) fprintf(prou,"\n");
            }
        }
        if ( intakt ) {
            printf(
            "\n\n  current scaling,  scf =  %11.4e clow = %12i eta =  %11.4e\n", 
            scf,clow,eta);
            
            if ( nres != 0 ) {
                printf(  "\n\n  scalres=\n");
                for (i = 1 ; i <= nres ; i++) {
                    printf(  "  %11.4e", w[i]);
                    if ( i % 6 == 0 || i == nres ) printf(  "\n");
                }
            }
        }
        return;
        
    case 13:
        if ( accinf[itstep][27] == zero ) {
            if ( intakt ) printf("update suppressed\n");
            fprintf(prou,        "update suppressed\n");
        } else if ( accinf[itstep][27] == -one ) {
            fprintf(prou,        "restart with scaled unit matrix\n");
            if ( intakt ) printf("restart with scaled unit matrix\n");
        } else {
            fprintf(prou,"BFGS-update\n");
            fprintf(prou," type = %14.6e\n", accinf[itstep][27]);
            fprintf(prou,"  ny  = %14.6e\n", accinf[itstep][28]);
            fprintf(prou," thet = %14.6e\n", accinf[itstep][29]);
            if ( intakt ) {
                printf(  "BFGS-update\n");
                printf(  " type = %14.6e\n", accinf[itstep][27]);
                printf(  "  ny  = %14.6e\n", accinf[itstep][28]);
                printf(  " thet = %14.6e\n", accinf[itstep][29]);
            }
        }
        return;
        
    case 14:
        if ( accinf[itstep][27] == zero ) {
            if ( intakt ) printf("update suppressed\n");
            fprintf(prou,        "update suppressed\n");
        } else if ( accinf[itstep][27] == one ) {
            fprintf(prou,"BFGS-update as in Pantoja and Mayne\n");
            fprintf(prou,"  tk  = %14.6e\n", accinf[itstep][28]);
            fprintf(prou," xsik = %14.6e\n", accinf[itstep][29]);
            if ( intakt ) {
                printf(  "BFGS-update\n");
                printf(  "  tk  = %14.6e\n", accinf[itstep][28]);
                printf(  " xsik = %14.6e\n", accinf[itstep][29]);
            }
        } else {
            if ( intakt ) printf("restart with scaled unit matrix\n");
            fprintf(prou,        "restart with scaled unit matrix\n");
        }
        return;
        
    case 15:
        if ( intakt ) printf("\n\n\n singular case : full regularized SQP\n");
        fprintf(prou,        "\n\n\n singular case : full regularized SQP\n");
        if ( intakt ) printf("  del = %.15e\n", del);
        fprintf(prou,        "  del = %.15e\n", del);
        
        if ( intakt ) { 
            printf(  "\n\n  scalres=\n");
            for (i = 1 ; i <= nres ; i++) {
                printf(  "  %11.4e", w[i]);
                if ( i % 6 == 0 || i == nres ) printf(  "\n");
            }
        }
        fprintf(prou,"\n\n  scalres=\n");
        for (i = 1 ; i <= nres ; i++) {
            fprintf(prou,"  %11.4e",w[i]);
            if ( i % 6 == 0 || i == nres ) fprintf(prou,"\n");
        }
        if ( te3 ) 
        {
            fprintf(prou,"gradients of binding constraints");
            for ( i =1 ; i <= alist[0] ; i++ )
            {
                l = alist[i];
                fprintf(prou,"\n\n gradient of restriction nr.%4i\n ", l);
                for (j = 0 ; j <= n ; j++)
                {
                    fprintf(prou," %11.4e  ", gres[j][l]);
                    if ( j % 5 == 0 || j == n ) fprintf(prou,"\n ");
                }
                
            }
        }
        return;
        
    case 16:
        fprintf(prou,"exit from full SQP\n");
        fprintf(prou,"            termination reason  %7.0e\n", 
        accinf[itstep][30]);
        fprintf(prou,"          final value of tauqp  %10.3e\n",
        accinf[itstep][31]);
        fprintf(prou,"      sum norm of slack vector  %10.3e\n",
        accinf[itstep][32]);
                                
        fprintf(prou,"\n\n  phase=%3i  scf0= %11.4e\n",phase,scf0);
        fprintf(prou,"  d = \n");
        for (i = 1 ; i <= n ; i++) {
            fprintf(prou,"  %11.4e", d[i]);
            if ( i % 6 == 0 || i == n ) fprintf(prou,"\n");
        }
        if ( nres != 0 ) { 
            fprintf(prou,"\n  multipliers: first estimate\n  u =\n");
            for (k = 1 ; k <= nres ; k++) {
                fprintf(prou," %4i  %11.4e", k,u[k]);       
                if ( k % 4 == 0 || k == nres ) fprintf(prou,"\n");
            }
        }
        if ( intakt ) {
            printf(  "exit from full SQP\n");
            printf(  "            termination reason  %7.0e\n", 
            accinf[itstep][30]);
            printf(  "          final value of tauqp  %10.3e\n",
            accinf[itstep][31]);
            printf(  "      sum norm of slack vector  %10.3e\n",
            accinf[itstep][32]);
                                
            printf("\n\n  phase=%3i  scf0= %11.4e\n",  phase,scf0);
            printf("  d = \n");
            for (i = 1 ; i <= n ; i++) {
                printf(  "  %11.4e", d[i]);
                if ( i % 6 == 0 || i == n ) printf(  "\n");
            }
            if ( nres != 0 ) {
                printf(  "\n  multipliers: first estimate\n  u =\n");
                for (k = 1 ; k <= nres ; k++) {
                    printf(  " %4i  %11.4e", k,u[k]);
                    if ( k % 4 == 0 || k == nres ) printf(  "\n");
                }
            }
        }
        return;
        
    case 17:
        fprintf(prou,"small directional derivative %.15e: finish\n",dirder);
        if ( intakt )
        printf(      "small directional derivative %.15e: finish\n",dirder);
        
        return;
        
    case 18:
        if ( intakt )
        printf(      "small correction from full regularized SQP,finish\n");
        fprintf(prou,"small correction from full regularized SQP,finish\n");
        
        return;
        
    case 19:
        fprintf(prou,        "QP-solver terminated unsuccessfully\n");
        if ( intakt ) printf("QP-solver terminated unsuccessfully\n");
        
        return;
        
    case 20:
        if ( intakt ) printf("restart with scaled unit matrix\n");
        fprintf(prou,        "restart with scaled unit matrix\n");
        
        return;
        
    case 21:
    
        return;
        
    case 22:
    
        return;
    }
}

/* **************************************************************************** */
/*         computation of new scaling factors for L1-penalty-function           */
/* **************************************************************************** */
void o8sce(void) {

    void    o8info(INTEGER icase);
    
    static INTEGER  i;
    static DOUBLE   term,s1,s2,diff0;
    static LOGICAL  wlow;

    wlow = FALSE;
    for (i = 1 ; i <= nres ; i++) {
    
        /* w1 tentative new weights */
        
        term = ny*fabs(u[i])+tau;
        if ( term > w[i] ) {
            w1[i] = term+tau;
        } else {
            w1[i] = w[i];
            if ( term < w[i]*p5 && bind[i] == 1 ) w1[i] = (term+w[i])*p5;
        }
        if ( w1[i] < w[i]) wlow = TRUE;
    }
    /* wlow equals TRUE if one tentative weight at least has been decreased */
    
    s1 = zero;
    s2 = zero;
    for (i = 1 ; i <= nres ; i++) {
        if ( i <= nh ) {
            s1 = s1+w1[i]*fabs(resst[i]);
            s2 = s2+w1[i]*fabs(res[i]);
        } else {
            s1 = s1-min(zero,resst[i])*w1[i];
            s2 = s2-min(zero,res[i]  )*w1[i];
        }
    }
    diff0 = (fxst-fx)*scf+(s1-s2);
    if ( wlow && diff0 >= eta*clow && itstep-lastdw > max(5,min(20,n/10)) ) {
    
        /* accept new (diminished ) weights */
        
        if ( clow > itstep/10 ) {
            eta = onep3*eta;
            if ( ! silent ) o8info(11);
        }
        lastch = itstep;
        lastdw = itstep;
        level  = diff0/iterma;
        psist  = s1;
        psi    = s2;
        for (i = 1 ; i <= nres ; i++) {
            w[i] = w1[i];
        }
        clow = clow+one;
    } else {
    
        /* increase individual weights if necessary. let weigths unchanged */
        /* otherwise                                                       */
        
        s1 = zero;
        s2 = zero;
        for (i = 1 ; i <= nres ; i++) {
            if ( w1[i] > w[i] ) {
                lastup = itstep;
                lastch = itstep;
            }
            w[i] = max(w[i],w1[i]);
            if ( i <= nh ) {
                s1 = s1+w[i]*fabs(resst[i]);
                s2 = s2+w[i]*fabs(res[i]);
            } else {
                s1 = s1-w[i]*min(zero,resst[i]);
                s2 = s2-w[i]*min(zero,res[i]);
            }
        }
        psist = s1;
        psi   = s2;
    }   
    term = zero;
    if ( nres >= 1 ) term = w[1];
    for (i = 2 ; i <= nres ; i++) {
        term = max(term,w[i]);
    }
    accinf[itstep][20] = term;
    
    /* maximum of weights */
    
    accinf[itstep][19] = clow;
    
    if ( ! silent ) o8info(12);
    
    return;
}

/* **************************************************************************** */
/*          computation of the Pantoja-Mayne BFGS-update of hessian             */
/* **************************************************************************** */
void o8bfgs(void) {

    void    o8msg (INTEGER num);
    void    o8inim(void);
    DOUBLE  o8sc1 (INTEGER i,INTEGER j,DOUBLE a[],DOUBLE b[]);
    DOUBLE  o8sc2 (INTEGER n,INTEGER m,INTEGER j,DOUBLE a[][NX+1],   DOUBLE b[]);
    DOUBLE  o8sc3 (INTEGER n,INTEGER m,INTEGER j,DOUBLE a[][NRESM+1],DOUBLE b[]);
    DOUBLE  o8sc3_(INTEGER n,INTEGER m,INTEGER j,DOUBLE a[][NX+1],   DOUBLE b[]);
    void    o8upd(DOUBLE r[][NX+1],DOUBLE z[],DOUBLE y[],INTEGER n,LOGICAL *fail);
    DOUBLE  o8vecn(INTEGER nl,INTEGER nm,DOUBLE x[]);
    
    static INTEGER  i,j;
    static DOUBLE   dg[NX+1],adx[NX+1],den1,den2,den3,
                    th,tk,xsik,ltdx[NX+1],gtdx[NRESM+1],updx[NX+1],updz[NX+1],
                    term,term1,anorm,acond,ndx,ngtdx,den21;
    static LOGICAL  fail;

    for ( i = 1 ;  i <= n ;  i++) {
    
        /* multiply dx = (s in the usual notation) by Cholesky-factor */
        /* stored in the upper half of a                              */
        
         ltdx[i] = o8sc2(i,n,i,a,difx);
         dg[i]   = gphi1[i]-gphi0[i];
    }
    if ( o8vecn(1,n,dg) == zero ) {
    
        /* suppress update */
        
        accinf[itstep][27] = zero;
        accinf[itstep][28] = zero;
        accinf[itstep][29] = zero;
        
        if ( ! silent ) o8msg(21);
        
        return;
    }
    for (i = 1 ; i <= n ; i++) {
        adx[i] = o8sc3_(1,i,i,a,ltdx);
    }
    /* adx = a * ( x-x0), x-x0 = difx */
    
    for (i = 1 ; i <= alist[0] ; i++) {
        gtdx[i] = o8sc3(1,n,alist[i],gres,difx);
        gtdx[i] = gtdx[i]/gresn[alist[i]];
    }
    /* gtdx = grad(res)(transp)*(x-x0) */
    
    ndx   = o8vecn(1,n,difx);
    tk    = min(p5,pow(dnorm,2));
    anorm = zero;
    term1 = fabs(a[1][1]);
    anorm = zero;
    for (i = 1 ; i <= n ; i++) {
        for (j = i ; j <= n ; j++) {
            anorm = anorm+pow(a[i][j],2);
        }
        term1 = min(term1,fabs(a[i][i]));
    }
    if ( term1 != zero ) {
        acond = anorm/pow(term1,2);
    } else {