donlp2.c

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

        }
    }
    for (i = nh+1 ; i <= nres ; i++) {
    
        /* modify del0, such that lower and upper bound never become binding */
        /* simultaneously                                                    */
        
        if ( gunit[1][i] == 1 ) {
            j = gunit[2][i];
            if ( og[j] < infiny && ug[j] > -infiny ) {
                del0 = min(del0,(og[j]-ug[j])*tm1*fabs(gunit[3][i]));
            }
        }
    }
    for (i = nh+1 ; i <= nres ; i++) {
    
        /* delfac corresponds to an indirect primal scaling */
        
        if ( gunit[1][i] == 1 ) {
            j = gunit[2][i];
            if ( gunit[3][i] > 0 ) {
                delfac[i] = max(delfac[i],fabs(ug[j])*tm1);
                if ( og[j] < infiny )
                delfac[i] = min(delfac[i],(og[j]-ug[j])/(tp1*del0));
            } else {
                delfac[i] = max(delfac[i],fabs(og[j])*tm1);
                if ( ug[j] > -infiny )
                delfac[i] = min(delfac[i],(og[j]-ug[j])/(tp1*del0));
            }
        }
    }
    bd0 = infiny;
    for (i = 1 ; i <= n ; i++) {
        if ( ug[i] > zero ) bd0 = min(bd0,og[i]);
        if ( og[i] < zero ) bd0 = min(bd0,-ug[i]);
    }
    /* change x if necessary , such that bounds not violated */
    
    corr = FALSE;
    
    /* evaluate gradients of supersimple functions only once              */
    /* a function is said to be supersimple if it is of the form a*x[j]+b */
    
    if ( gunit[1][0] == 1 ) {
        gconst[0] = TRUE;
        val[0]    = TRUE;
        for (i = 1 ; i <= n ; i++) {
            gradf[i] = zero;
        }
        gradf[gunit[2][0]] = gunit[3][0]*xsc[gunit[2][0]];
        gfn = fabs(gunit[3][0]);
    } else {
        val[0] = FALSE;
        for (i = 1 ; i <= n ; i++) {
            gradf[i] = zero;
        }
    }
    for (i = 1 ; i <= nh ; i++) {
        if ( gunit[1][i] == 1 ) {
        
            /* a fixed variable. corrected if necessary */
            
            val[i]    = TRUE;
            gconst[i] = TRUE;
            gres[gunit[2][i]][i] = gunit[3][i]*xsc[gunit[2][i]];
            gresn[i] = fabs(gunit[3][i])*xsc[gunit[2][i]];
            if ( gresn[i] == zero ) {
                if ( ! silent ) 
                fprintf(meu,"%i fixed variable/zero gradient\n", gunit[2][i]);
                fclose(meu);
                fclose(prou);
                exit(1);
            }
            esh(i,xnull,&hxi);
            
            term = -hxi/gunit[3][i];
            if ( term != x[gunit[2][i]] ) corr = TRUE;
            x[gunit[2][i]] = term;
        }
    }
    if ( bloc ) user_eval(x,0);
    
    for (i = nh+1 ; i <= nres ; i++) {
        if ( gunit[1][i] == 1 ) {
            if ( gunit[3][i] == 0 ) {
                if ( ! silent ) 
                fprintf(meu,"%i bounded variable, zero gradient\n",gunit[2][i]);
                fclose(meu);
                fclose(prou);
                exit(1);
            }
            esg(i-nh,x,&gxi);
            
            gxi = two*delmin-gxi;
            if ( gxi > zero ) {
                corr = TRUE;
                x[gunit[2][i]] = x[gunit[2][i]]+gxi/gunit[3][i];
            }
            gres[gunit[2][i]][i] = gunit[3][i]*xsc[gunit[2][i]];
            gresn[i]  = fabs(gunit[3][i])*xsc[gunit[2][i]];
            val[i]    = TRUE;
            gconst[i] = TRUE;
        }
    }
    if ( corr && ! silent ) o8msg(13);

    /* remember to reevaluate the functions if corr = TRUE and bloc = TRUE */

    if ( bloc ) user_eval(x,1);

    for (i = 1 ; i <= nres ; i ++) {
        bind[i]  = 0;
        bind0[i] = 0;
        u[i]     = zero;
        u0[i]    = zero;
        cres[i]  = 0;
        cgres[i] = 0;
        
        /* initial weights of the penalty-term */
        
        if ( cold ) w[i] = one;
        sort[i] = i;
    }   
    clow = one;
    ny   = two;

    /* scf  = weight factor for objective function    */
    /* scf0 = damping factor for tangential direction */
    
    scf    = one;
    scf0   = one;
    sigla  = twop11;
    beta   = four;  /* formerly two */
    alpha  = tm1;
    delta1 = p9;
    delta  = tm3;   /* delta = tm2  formerly */
    theta  = p9;    /* theta = 0.99 formerly */
    icf    = 0;
    icgf   = 0;
    if ( ! silent ) {
        fprintf(prou,"donlp2, v3, 05/29/98, copyright P. Spellucci\n");
        
        time(&tim);
        
        fprintf(prou,"%s",ctime(&tim));
        fprintf(prou,"%s\n", name);
        
        fprintf(meu, "donlp2, v3, 05/29/98, copyright P. Spellucci\n");
        fprintf(meu, "%s",ctime(&tim));
        fprintf(meu, "%s\n", name);
    }
    return;
}

/* **************************************************************************** */
/*                        do final solution check and output                    */
/* **************************************************************************** */
void o8fin(void) {

    void    solchk(void);

    static INTEGER  i,j,k,ih1,ih2,ih3,ih5,ih6,ih7,ih8,ih9;
    static DOUBLE   umin,term;
    static INTEGER  nsing,crtot,cgrtot,nupreg,nbfgs1,nresta;
    static char     line[65];

    /* termination reason + 11 = message number */
    
    static char *messag[] = {
    "",     /* not used : index 0 */
    "constraint evaluation returns error with current point",
    "objective evaluation returns error with current point",
    "QPsolver: working set singular in dual extended QP ",
    "QPsolver: extended QP-problem seemingly infeasible ",
    "QPsolver: no descent direction from QP for tau=tau_max",
    "QPsolver: on exit correction small, infeasible point",
    "stepsizeselection: computed d not a direction of descent",
    "more than MAXIT iteration steps",
    "stepsizeselection: no acceptable stepsize in [sigsm,sigla]",
    "stepsizeselection: directional deriv. very small, infeasible",
    "KT-conditions satisfied, no further correction computed",
    "KT-conditions satisfied, computed correction small",
    "stepsizeselection: x (almost) feasible, dir. deriv. very small",
    "KT-conditions (relaxed) satisfied, singular point",
    "very slow primal progress, singular or illconditoned problem",
    "very slow progress in x, singular problem",
    "correction very small, almost feasible but singular point",
    "max(n,10) small differences in penalty function,terminate"
    };

    if ( scf != zero ) {
        for (i = 1 ; i <= nres ; i++) {
            u[i] = u[i]/scf;
        }
    }
    /* in solchk the user might add some additional checks and/or output */
    
    solchk();

    if ( silent && ! intakt ) return;

    if ( intakt && ! silent ) printf("%s\n", name);
    
    if ( ! silent ) {
        if ( intakt ) {
            printf(  "\n     n= %9i    nh= %9i    ng= %9i\n", n,nh,ng);
            printf(  "\n  epsx= %9.3e sigsm= %9.3e\n"       , epsx,sigsm);
            printf(  "\nstartvalue\n");
            for (i = 1 ; i <= n ; i++) {
                printf(  " %14.7e ", xst[i]);
                if ( i % 5 == 0 || i == n ) printf(  "\n");
            }
        }
        fprintf(prou,"\n     n= %9i    nh= %9i    ng= %9i\n", n,nh,ng);
        fprintf(prou,"\n  epsx= %9.3e sigsm= %9.3e\n"       , epsx,sigsm);
        fprintf(prou,"\nstartvalue\n");
        for (i = 1 ; i <= n ; i++) {
            fprintf(prou," %14.7e ", xst[i]);
            if ( i % 5 == 0 || i == n ) fprintf(prou,"\n");
        }
    }
    if ( intakt && ! silent ) {
        printf("\n  eps= %9.2e  tol= %9.2e del0= %9.2e delm= %9.2e tau0= %9.2e\n",
        epsmac,tolmac,del0,delmin,tau0);
        printf(  "  tau= %9.2e   sd= %9.2e   sw= %9.2e  rho= %9.2e rho1= %9.2e\n",
        tau,smalld,smallw,rho,rho1);
    }
    if ( ! silent ) { 
        fprintf(prou,
               "\n  eps= %9.2e  tol= %9.2e del0= %9.2e delm= %9.2e tau0= %9.2e\n",
        epsmac,tolmac,del0,delmin,tau0);
        fprintf(prou,
                 "  tau= %9.2e   sd= %9.2e   sw= %9.2e  rho= %9.2e rho1= %9.2e\n",
        tau,smalld,smallw,rho,rho1);
    }
    if ( ! silent ) {
        fprintf(prou,
        " scfm= %9.2e  c1d= %9.2e epdi= %9.2e\n  nre= %9i anal= %9i\n",
        scfmax,c1d,epsdif,nreset,analyt);
        if ( intakt ) printf(
        " scfm= %9.2e  c1d= %9.2e epdi= %9.2e\n  nre= %9i anal= %9i\n",
        scfmax,c1d,epsdif,nreset,analyt);
    }
    if ( ! silent && ! analyt ) {
        fprintf(prou," vbnd= %9.2e efcn= %9.2e diff=%1i\n"
        , taubnd,epsfcn,difftype);
        if ( intakt ) printf("taubnd= %9.2e epsfcn= %9.2e difftype=%1i\n"
        , taubnd,epsfcn,difftype);
    }
    i    = 0;
    j    = 0;
    umin = zero;
    for (k = 1 ; k <= nres ; k++) {
        i = i+ cres[k];
        j = j+cgres[k];
        if ( k > nh ) umin = min(umin,u[k]);
    }
    crtot  = i;
    cgrtot = j;
    nsing  = 0;
    nresta = 0;
    nupreg = 0;
    nbfgs1 = 0;
    for (k = 1 ; k <= itstep ; k++) {
        if ( accinf[k][10] == one ) nsing  = nsing+1;
        if ( accinf[k][27] == one ) nbfgs1 = nbfgs1+1;
        
        /* for the Pantoja Mayne update */
        
        if ( accinf[k][29] == zero && accinf[k][27] == one ) nupreg = nupreg+1;
        if ( accinf[k][27] == -one ) nresta = nresta+1;
    }
    k = (int)optite+11;
    if ( k >= 1 && k <= 18 ) {
        strcpy(line,messag[k]);
    } else {
        strcpy(line,"variable optite undefined on exit");
    }
    if ( intakt && ! silent ) {
        printf(      "\n termination reason:\n %s\n",line);
        printf(        " evaluations of f                    %9i\n",    icf);
        printf(        " evaluations of grad f               %9i\n",    icgf);
        printf(        " evaluations of constraints          %9i\n",    crtot);
        printf(        " evaluations of grads of constraints %9i\n",    cgrtot);
        printf(        " final scaling of objective          %13.6e\n", scf);
        printf(        " norm of grad(f)                     %13.6e\n", gfn);
        printf(        " lagrangian violation                %13.6e\n", b2n);
        printf(        " feasibility violation               %13.6e\n", upsi);
        printf(        " dual feasibility violation          %13.6e\n", umin);
        printf(        " optimizer runtime sec's             %13.6e\n", runtim);
    }
    if ( ! silent ) {
        fprintf(prou,"\n termination reason:\n %s\n",line);
        fprintf(prou,  " evaluations of f                    %9i\n",    icf);
        fprintf(prou,  " evaluations of grad f               %9i\n",    icgf);
        fprintf(prou,  " evaluations of constraints          %9i\n",    crtot);
        fprintf(prou,  " evaluations of grads of constraints %9i\n",    cgrtot);
        fprintf(prou,  " final scaling of objective          %13.6e\n", scf);
        fprintf(prou,  " norm of grad(f)                     %13.6e\n", gfn);
        fprintf(prou,  " lagrangian violation                %13.6e\n", b2n);
        fprintf(prou,  " feasibility violation               %13.6e\n", upsi);
        fprintf(prou,  " dual feasibility violation          %13.6e\n", umin);
        fprintf(prou,  " optimizer runtime sec's             %13.6e\n", runtim);
    }                                         
    if ( intakt && ! silent ) printf("\n\n optimal value of f =  %21.14e\n", fx);
    if ( ! silent ) fprintf(prou,    "\n\n optimal value of f =  %21.14e\n", fx);
    
    if ( intakt && ! silent ) {
        printf(      "\n optimal solution  x =\n");
        for (i = 1 ; i <= n ; i++) {
            printf(      " %21.14e", x[i]);
            if ( i % 3 == 0 || i == n ) printf(     "\n");
        }
    }
    if ( ! silent ) {
        fprintf(prou,"\n optimal solution  x =\n");
        for (i = 1 ; i <= n ; i++) {
            fprintf(prou," %21.14e", x[i]);
            if ( i % 3 == 0 || i == n ) fprintf(prou,"\n");
        }
    }
    if ( nres != 0 && ! silent ) {
        fprintf(prou,"\n  multipliers are relative to scf=1\n");
        fprintf(prou,  "  nr.    constraint      normgrad (or 1)   multiplier\n");
        for (i = 1 ; i <= nres ; i++) {
            fprintf(prou," %4i  %14.7e   %14.7e   %14.7e\n"
            , i,res[i],gresn[i],u[i]);
        }
        if ( intakt ) {
            printf(  "\n  multipliers are relative to scf=1\n");
            printf(    "  nr.    constraint      normgrad (or 1)   multiplier\n");
            for (i = 1 ; i <= nres ; i++) {
                printf(  " %4i  %14.7e   %14.7e   %14.7e\n"
                , i,res[i],gresn[i],u[i]);
            }
        }
    }
    if ( nres != 0 && ! silent ) {
        fprintf(prou,"\n evaluations of restrictions and their gradients\n");
        for (i = 1 ; i <= nres ; i++) {
            fprintf(prou," (%6i,%6i)", cres[i],cgres[i]);
            if ( i % 5 == 0 || i == nres ) fprintf(prou,"\n");
        }
        if ( intakt ) {
            printf(  "\n evaluations of restrictions and their gradients\n");
            for (i = 1 ; i <= nres ; i++) {
                printf(  " (%6i,%6i)", cres[i],cgres[i]);
                if ( i % 5 == 0 || i == nres ) printf(  "\n");
            }
        }
    }
    if ( itstep > 1 && optite == 0 ) itstep = itstep-1;
    if ( ! silent )
    fprintf(prou,    "last estimate of condition of active gradients %10.3e\n",
    accinf[itstep][13]);
                    
    term = accinf[itstep][14];
    i    = itstep;
    while ( i > 1 && term == zero ) {
        i    = i-1;
        term = accinf[i][14];
    }
    if ( ! silent ) {
        fprintf(prou,"last estimate of condition of approx. hessian  %10.3e\n",