donlp2.c

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

        term);
        fprintf(prou,"iterative steps total           %5i\n", itstep);
        fprintf(prou,"# of restarts                   %5i\n", nresta);
        fprintf(prou,"# of full regular updates       %5i\n", nupreg);
        fprintf(prou,"# of updates                    %5i\n", nbfgs1);
        fprintf(prou,"# of full regularized SQP-steps %5i\n", nsing);

        if ( intakt ) {
            printf(  "last estimate of cond.nr. of active gradients  %10.3e\n",
            accinf[itstep][13]);
            printf(  "last estimate of cond.nr. of approx.  hessian  %10.3e\n",
            accinf[itstep][14]);
            printf(  "iterative steps total           %5i\n", itstep);
            printf(  "# of restarts                   %5i\n", nresta);
            printf(  "# of full regular updates       %5i\n", nupreg);
            printf(  "# of updates                    %5i\n", nbfgs1);
            printf(  "# of regularized full SQP-steps %5i\n", nsing);
        }
    }
    if ( optite < zero ) te1 = TRUE;
    if ( silent )        te1 = FALSE;
    if ( te1 ) {
        for (i = 1 ; i <= itstep ; i++) {
            ih1 = accinf[i][1];
            ih2 = accinf[i][9];
            ih3 = accinf[i][10];
            ih5 = accinf[i][18];
            ih6 = accinf[i][19];
            ih7 = accinf[i][22];
            ih8 = accinf[i][26];
            ih9 = accinf[i][27];
            
            fprintf(prou, 
             "%4i  fx= %13.6e scf= %13.6e psi= %13.6e ups= %13.6e\n",
            ih1,accinf[i][2],accinf[i][3],accinf[i][4],accinf[i][5]);
            fprintf(prou,
            "     del= %13.6e b20= %13.6e b2n= %13.6e  nr=%5i\n",
            accinf[i][6],accinf[i][7],accinf[i][8],ih2);
            fprintf(prou,
            "      si=%4i            u-= %13.6e c-r= %13.6e c-d= %13.6e\n",
            ih3,accinf[i][11],accinf[i][13],accinf[i][14]);
            fprintf(prou,
            "      xn= %13.6e  dn= %13.6e pha=%4i            cl=%14i\n",
            accinf[i][16],accinf[i][17],ih5,ih6);
            fprintf(prou,
            "     skm= %13.6e sig= %13.6e cf+=%5i          dir= %13.6e\n",
            accinf[i][20],accinf[i][21],ih7,accinf[i][23]);
            fprintf(prou,
            "     dsc= %13.6e cos= %13.6e vio=%5i\n",
            accinf[i][24],accinf[i][25],ih8);
            fprintf(prou,
            "     upd=%5i           tk= %13.6e xsi= %13.6e\n",
            ih9,accinf[i][28],accinf[i][29]);
            
            if ( accinf[i][10] == 1. ) {
                fprintf(prou,"     qpt= %13.0e tqp= %13.6e sla= %13.6e\n",
                accinf[i][30],accinf[i][31],accinf[i][32]);
            }
        }
    }
    /*  accinf a c c u m u l a t e d   i n f o r m a t i o n                    */
    /*  on iteration sequence                                                   */
    /*  1: step-nr                                                              */
    /*  2: f(x-k) current value of objective (zero in feasibility improvement   */
    /*            phase (-1) )                                                  */
    /*  3: scf    internal scaling of objective (zero in phase -1)              */
    /*  4: psi    the weighted penalty-term                                     */
    /*  5: upsi   the unweighted penalty-term (L1-norm of constraint vector)    */
    /*  6: del_k_1 bound for currently active constraints                       */
    /*  7: b2n0   L2-norm of projected gradient, based on constraints in level  */
    /*            delmin and below, measured in the norm induced by the         */
    /*            inverse hessian estimate                                      */
    /*  8: b2n    L2-norm of projected gradient based on del_k_1                */
    /*  9: nr     number of binding constraints                                 */
    /* 10: sing   if 1, the binding constraints don't satisfy the regularity    */
    /*            condition                                                     */
    /* 11: umin   infinity norm of negative part of multiplier                  */
    /* 12: -------------                                                        */
    /* 13: cond_r condition number of diagonal part of QR-decomp. of normalized */
    /*            gradients of binding constraints                              */
    /* 14: cond_h condition number of diagonal of Cholesky-factor of updated    */
    /*            full hessian                                                  */
    /* 15: scf0   the relative damping of tangential component if upsi>tau0/2   */
    /* 16: xnorm  L2-norm of x                                                  */
    /* 17: dnorm  L2-norm of d (correction from QP -subproblem, unscaled)       */
    /* 18: phase  -1 : infeasibility improvement phase, 0: initial optimization */
    /*            1  : binding constraints unchanged , 2: d small, Maratos      */
    /*                 correction in use                                        */
    /* 19: c_k    number of decreases of penalty weights                        */
    /* 20: wmax   infinity norm of weights                                      */
    /* 21: sig_k  stepsize from unidimensional minimization (backtracking)      */
    /* 22: cfincr number of objective evaluations for stepsize-algorithm        */
    /* 23: dirder directional derivative of penalty-function along d (scaled)   */
    /* 24: dscal  scaling factor for d                                          */
    /* 25: cosphi cos of arc between d and d_previous. if larger theta , sig    */
    /*            larger than one (up to sigla) is tried                        */
    /* 26: violis[0] number of constraints not binding at x but hit during      */
    /*            line search                                                   */
    /* 27:        type of update for h: 1 normal P&M-BFGS-update,               */
    /*            0 update suppressed, -1 restart with scaled unit matrix ,     */
    /*            2 standard BFGS, 3 BFGS modified by Powells device            */
    /* 28: ny_k/tk modification factor for damping the projector in BFGS/       */
    /*            Pantoja-Mayne-term respectively                               */
    /* 29: 1-my_k/xsik modification factor for damping the quasi-Newton-        */
    /*            relation in BFGS for unmodified BFGS ny_k should be larger    */
    /*            than updmy0 (near one) and 1-my_k equal one./Pantoja-Mayne    */
    /*            term respectively                                             */
    /* 30: qpterm 0, if sing = -1, termination indicator of QP-solver otherwise */
    /*            1: successful, -1: tau becomes larger than tauqp without      */
    /*            slack-variables becoming sufficiently small                   */
    /*            -3: working set of QP-solver becomes linearly dependent       */
    /*            -2: infeasible QP-problem (theoretically impossible)          */
    /* 31: tauqp  weight  of slack-variables in QP-solver                       */
    /* 32: infeas L1-norm of slack-variables in QP-solver                       */

    if ( ! silent ) fclose(prou);
    if ( ! silent ) fclose(meu);
    
    return;
}
/* **************************************************************************** */
/*                      prints informations if te2 = TRUE                       */
/* **************************************************************************** */
void o8info(INTEGER icase) {

    void o8mdru (DOUBLE a[][NRESM+1],INTEGER ma,INTEGER na,char head[],
                 FILE *lognum,LOGICAL fix);
    void o8mdru_(DOUBLE a[][NX+1],   INTEGER ma,INTEGER na,char head[],
                 FILE *lognum,LOGICAL fix);
            
    static INTEGER  i,j,l,k;
    static DOUBLE   y,phih;
    static char     head[41];

    if(!te2) return;
    
    switch (icase) {
        
    case 1:
        fprintf(prou,"\n\n\n");
        for (i = 1 ; i <= 80 ; i++) fprintf(prou,"=");
        fprintf(prou,"\n          %4i-th iteration step\n", itstep);
        fprintf(prou,"   scf= %11.4e psist= %11.4e   psi= %11.4e  upsi= %11.4e\n", 
        scf,psist,psi,upsi);
        fprintf(prou,"  fxst= %11.4e    fx= %11.4e\n", fxst,fx);
        fprintf(prou,"  x=\n");
        for (i = 1 ; i <= n ; i++) {
            fprintf(prou,"  %11.4e", x[i]);
            if ( i % 6 == 0 || i == n ) fprintf(prou,"\n");
        }
        fprintf(prou," valid permutation of x\n\n");
        
        for (i = 1 ; i <= n ; i++) {
            fprintf(prou,"%3i ", perm[i]);
            if ( i % 20 == 0 || i == n ) fprintf(prou,"\n");
        }
        if ( phase >= 0 && te3 ) {
            strcpy(head,"quasi-Newton-matrix full update");
            
            o8mdru_(a,n,n,head,prou,FALSE);
        }
        if ( intakt ) {
            printf(  "\n\n\n");
            for (i = 1 ; i <= 80 ; i++) printf(  "=");
            printf(  "\n          %4i-th iteration step\n", itstep);
            printf(  "   scf= %11.4e psist= %11.4e   psi= %11.4e  upsi= %11.4e\n", 
            scf,psist,psi,upsi);
            printf(  "  fxst= %11.4e    fx= %11.4e\n", fxst,fx);
            printf(  "  x=\n");
            for (i = 1 ; i <= n ; i++) {
                printf("  %11.4e", x[i]);
                if ( i % 6 == 0 || i == n ) printf(  "\n");
            }
            printf(  " valid permutation of x\n\n");
        
            for (i = 1 ; i <= n ; i++) {
                printf(  "%3i ", perm[i]);
                if ( i % 20 == 0 || i == n ) printf(  "\n");
            }
        }
        return;
        
    case 2:
      fprintf(prou,"\n\n  del= %12.5e  b2n0= %12.5e   b2n= %12.5e   gfn= %12.5e\n", 
      del,b2n0,b2n,gfn);
        
        if ( alist[0] != 0 ) {
            fprintf(prou,"\n\n values of restrictions\n ");
            for (i = 1 ; i <= alist[0] ; i++) {
                fprintf(prou,"(%4i   %11.4e   %11.4e)  ", 
                alist[i],res[alist[i]],gresn[alist[i]]);
                if ( i % 2 == 0 || i == alist[0] ) fprintf(prou,"\n ");
            }
        }
        if ( alist[0] != 0 && ! singul ) {
            fprintf(prou,"\n\n   diag[r]=\n");
            for (i = 1 ; i <= alist[0] ; i++) {
                fprintf(prou,"  %11.4e", diag[i]);
                if ( i % 6 == 0 || i == alist[0] ) fprintf(prou,"\n");
            }
        }
        if ( alist[0] != 0 && te3 ) {
            for (i = 1 ; i <= alist[0] ; i++) {
                l = alist[i];
                fprintf(prou,"\n\n gradient of restriction nr.%4i\n ", l);
                for (j = 1 ; j <= n ; j++) {
                    fprintf(prou," %11.4e  ", gres[j][l]);
                    if ( j % 5 == 0 || j == n ) fprintf(prou,"\n ");
                }
            }
        }
        if ( intakt ) {
            printf("\n\n  del= %12.5e  b2n0= %12.5e   b2n= %12.5e   gfn= %12.5e\n", 
            del,b2n0,b2n,gfn);
        
            if ( alist[0] != 0 ) {
                printf(  "\n\n values of restrictions\n ");
                for (i = 1 ; i <= alist[0] ; i++) {
                    printf(  "(%4i   %11.4e   %11.4e)  ", 
                    alist[i],res[alist[i]],gresn[alist[i]]);
                    if ( i % 2 == 0 || i == alist[0] ) printf(  "\n ");
                }
            }
            if ( alist[0] != 0 && ! singul ) {
                printf(  "\n\n   diag[r]=\n");
                for (i = 1 ; i <= alist[0] ; i++) {
                    printf(  "  %11.4e", diag[i]);
                    if ( i % 6 == 0 || i == alist[0] ) printf(  "\n");
                }
            }
            if ( alist[0] != 0 && te3 ) {
                for (i = 1 ; i <= alist[0] ; i++) {
                    l = alist[i];
                    printf(  "\n\n gradient of restriction nr.%4i\n ", l);
                    for (j = 1 ; j <= n ; j++) {
                        printf(  " %11.4e  ", gres[j][l]);
                        if ( j % 5 == 0 || j == n ) printf(  "\n ");
                    }
                }
            }
        }
        return;
        
    case 3:
        if( ! (nr == 0 || phase == -1) ) {
            fprintf(prou,"\n  multipliers: first estimate\n  u =\n");
            for (k = 1 ; k <= nr ; k++) {
                fprintf(prou," %4i  %11.4e", alist[k],u[alist[k]]);
                if ( k % 4 == 0 || k == nr ) fprintf(prou,"\n");
            }
        }
        if( ! (nr == 0 || phase == -1) && intakt ) {
            printf(      "\n  multipliers: first estimate\n  u =\n");
            for (k = 1 ; k <= nr ; k++) {
                printf(      " %4i  %11.4e", alist[k],u[alist[k]]);
                if ( k % 4 == 0 || k == nr ) printf(      "\n");
            }
        }
        return;
        
    case 4:
        if( ! (nr == 0 || phase == -1) ) {
            fprintf(prou,"\n  multipliers: second estimate\n  u =\n");
            for (k = 1 ; k <= nr ; k++) {
                fprintf(prou," %4i  %11.4e", alist[k],u[alist[k]]);
                if ( k % 4 == 0 || k == nr ) fprintf(prou,"\n");
            }
        }
        if( ! (nr == 0 || phase == -1) && intakt ) {
            printf(      "\n  multipliers: second estimate\n  u =\n");
            for (k = 1 ; k <= nr ; k++) {
                printf(      " %4i  %11.4e", alist[k],u[alist[k]]);
                if ( k % 4 == 0 || k == nr ) printf(      "\n");
            }
        }
        return;
        
    case 5:
        if( intakt )
        printf(          "  condition number of r     %.15e\n",accinf[itstep][13]);
        fprintf(prou,    "  condition number of r     %.15e\n",accinf[itstep][13]);
        if ( phase == -1 ) {
        
            return;
            
        } else {
            fprintf(prou,"  condition number of a     %.15e\n",accinf[itstep][14]);
            if ( intakt ) {
                printf(  "  condition number of a     %.15e\n",accinf[itstep][14]);
            }
            return;
        }
    case 6:
    
        return;
        
    case 7:
        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 ( phase == 2 ) {
            fprintf(prou,"\n\n  dd=\n");
            for (i = 1 ; i <= n ; i++) {
                fprintf(prou,"  %11.4e", dd[i]);
                if ( i % 6 == 0 || i == n ) fprintf(prou,"\n");
            }
        }
        if ( intakt ) {
            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 ( phase == 2 ) {
                printf(  "\n\n  dd=\n");
                for (i = 1 ; i <= n ; i++) {
                    printf(  "  %11.4e", dd[i]);
                    if ( i % 6 == 0 || i == n ) printf(  "\n");
                }
            }
        }
        return;
        
    case 8:
        y    = tau0*p5;
        phih = fx*scf+psi;
    
        if ( intakt ) {
            printf(  "\n\n start unimin\n\n");
            printf(  "    phi= %11.4e   dphi= %11.4e    psi= %11.4e tau0/2= %11.4e\n",
            phih,dirder,psi,y);
            printf(  "     fx= %11.4e  dscal= %11.4e    scf= %11.4e   upsi= %11.4e\n", 
            fx,dscal,scf,upsi);
        }
        fprintf(prou,"\n\n start unimin\n\n");
        fprintf(prou,"    phi= %11.4e   dphi= %11.4e    psi= %11.4e tau0/2= %11.4e\n",
        phih,dirder,psi,y);
        fprintf(prou,"     fx= %11.4e  dscal= %11.4e    scf= %11.4e   upsi= %11.4e\n",
        fx,dscal,scf,upsi);
        
        return;
        
    case 9:
        fprintf(prou,"    sig= %11.4e     fx= %11.4e    psi= %11.4e   upsi= %11.4e\n",
        sig,fx1,psi1,upsi1);
        if ( intakt ) 
        printf(      "    sig= %11.4e     fx= %11.4e    psi= %11.4e   upsi= %11.4e\n", 
        sig,fx1,psi1,upsi1);
        
        return;
        
    case 10:
        fprintf(prou,"\n\n end unimin\n");
        fprintf(prou,  "\n sig= %11.4e  num. f-evaluations%2i\n", sig,cfincr);
        fprintf(prou,    " list of inactive hit constraints\n");
        for (i = 1 ; i <= violis[0] ; i++) {
            fprintf(prou,"%4i  ", violis[i]);
            if ( i % 13 == 0 || i == violis[0] ) fprintf(prou,"\n");
        }
        if ( violis[0] == 0 ) fprintf(prou,"none\n");