amplqp.c

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

/* **************************************************************************** */
/*                                 user functions                               */
/* **************************************************************************** */
#include "o8para.h"

main() {
    void donlp2(void);
    
    donlp2();
    
    exit(0);
}

/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    static INTEGER  i;
    
    /* name is ident of the example/user and can be set at users will       */
    /* the first static character must be alphabetic. 40 characters maximum */

    strcpy(name,"amplexam");
     
    n  = 5;
    nh = 2;
    ng = 1;
    for (i = 1 ; i <= 5 ; i++) {
        x[i] = 0.e0;
    }
    for (i = 0 ; i <= 3 ; i++) {
        gunit[1][i] = -1;
        gunit[2][i] = 0;
        gunit[3][i] = 0;
    }
    gconst[0] = FALSE;
    for (i = 1 ; i <= 3 ; i++) {
        gconst[i] = TRUE;
    }
    tau0 = 1.e8;
    del0 = 1.e0;
    tau  = 0.1e0;

    /* since we use high order numerical differentiation : */
    
    analyt = TRUE;
    epsdif = 0.e0;
    cold   = TRUE;
    silent = FALSE;

    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    epsx = 1.e-8;

    return;
}

/* **************************************************************************** */
/*  the user may add additional computations using the computed solution here   */
/* **************************************************************************** */
void solchk(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    #include "o8cons.h"

    return;
}

/* **************************************************************************** */
/*                               objective function                             */
/* **************************************************************************** */
void ef(DOUBLE x[],DOUBLE *fx) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X

    icf = icf+1;
    *fx = 
         0.5*(0.371136*x[1]*x[1]) + 0.5*(5.14041*x[1]*x[2])
         + 0.5*(2.60871*x[1]*
         x[3]) + 0.5*(-1.13012*x[1]*x[4])
         + 0.5*(-6.05446*x[1]*x[5]) + 0.5*(
         5.14041*x[2]*x[1]) + 0.5*(72.1133*x[2]*x[2])
         + 0.5*(36.1319*x[2]*x[3])
         + 0.5*(-15.6528*x[2]*x[4])
         + 0.5*(-83.8572*x[2]*x[5]) + 0.5*(2.60871*
         x[3]*x[1]) + 0.5*(36.1319*x[3]*x[2]) 
         + 0.5*(18.4058*x[3]*x[3]) + 0.5*( -7.94363*x[3]*x[4])
         + 0.5*(-42.5567*x[3]*x[5]) + 0.5*(-1.13012*x[4]*
         x[1]) + 0.5*(-15.6528*x[4]*x[2]) + 0.5*(-7.94363*x[4]*x[3])
         + 0.5*( 4.32868*x[4]*x[4]) + 0.5*(23.8574*x[4]*x[5])
         + 0.5*(-6.05446*x[5]*x[1])
         + 0.5*(-83.8572*x[5]*x[2]) + 0.5*(-42.5567*x[5]*x[3])
         + 0.5*(23.8574*  x[5]*x[4]) + 0.5*(132.59*x[5]*x[5])
         + 8.61104*x[1] + 14.8321*x[2] +
           3.92921*x[3] + 1.43377*x[4] + 5.98614*x[5];
    
    return;
}

/* **************************************************************************** */
/*                          gradient of objective function                      */
/* **************************************************************************** */
void egradf(DOUBLE xl[],DOUBLE gradxl[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X
    #include "o8cons.h"
    
    void ef(DOUBLE x[],DOUBLE *fx);
    
    static DOUBLE   fxl[7];
    static DOUBLE   sd1,sd2,sd3,xincr,zz,d1,d2,d3;
    static INTEGER  i;
    
    /* ************************************************* */
    /* high precision numerical differentiation          */
    /* by sixth order extrapolation                      */
    /* uses Richardson-extrapolation of three values     */
    /* of the symmetric difference quotient with         */
    /* rather large discretization stepsize              */
    /* the optimal stepsize depends on the seventh       */
    /* partial derivatives of the function . since these */
    /* are inaccessible, we replace them by one in the   */
    /* estimate of deldif                                */
    /* ************************************************* */
    
    deldif = exp(log(epsmac)/seven)*twom2;
    icgf   = icgf+1;
    for (i = 1 ; i <= n ; i++) {
        zz    = xl[i];
        xincr = deldif*(fabs(xl[i])+one);
        xl[i] = zz-xincr;
        
        ef(xl,&fxl[1]);
        
        xl[i] = zz+xincr;
        
        ef(xl,&fxl[2]);
        
        xincr = xincr+xincr;
        d1    = xincr;
        xl[i] = zz-xincr;
        
        ef(xl,&fxl[3]);
        
        xl[i] = zz+xincr;
        
        ef(xl,&fxl[4]);
        
        xincr = xincr+xincr;
        d2    = xincr;
        xl[i] = zz-xincr;
        
        ef(xl,&fxl[5]);
        
        xl[i] = zz+xincr;
        
        ef(xl,&fxl[6]);
        
        xl[i]     = zz;
        d3        = xincr+xincr;
        sd1       = (fxl[2]-fxl[1])/d1;
        sd2       = (fxl[4]-fxl[3])/d2;
        sd3       = (fxl[6]-fxl[5])/d3;
        sd3       = sd2-sd3;
        sd2       = sd1-sd2;
        sd3       = sd2-sd3;
        gradxl[i] = sd1+p4*sd2+sd3/c45;
    }
    return;
}

/* **************************************************************************** */
/*                compute the i-th equality constaint, value is hxi             */
/* **************************************************************************** */
void eh(INTEGER i,DOUBLE x[],DOUBLE *hxi) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X
    
    cres[i] = cres[i]+1;
    
    if ( i == 1 ) {
    
        /* s.t. eq(1): */
        
        *hxi = 
         0.302736*x[1] + 0.217135*x[2] + 2.66167*x[3] - 1.54133*x[4]
         -  1.12413*x[5] -( -3.50618);
    } else {
        
        /* s.t. eq(2): */

        *hxi = -2.035*x[1] - 1.23544*x[4] -( -2.39339);
    }
    return;
}

/* **************************************************************************** */
/*              compute the gradient of the i-th equality constraint            */
/* **************************************************************************** */
void egradh(INTEGER j,DOUBLE xl[],DOUBLE gradxl[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X
    #include "o8cons.h"
    
    void eh(INTEGER i,DOUBLE x[],DOUBLE *hxi);
    
    static DOUBLE   fxl[7];
    static DOUBLE   sd1,sd2,sd3,xincr,zz,d1,d2,d3;
    static INTEGER  i;
    
    /* ************************************************* */
    /* high precision numerical differentiation          */
    /* by sixth order extrapolation                      */
    /* uses Richardson-extrapolation of three values     */
    /* of the symmetric difference quotient with         */
    /* rather large discretization stepsize              */
    /* the optimal stepsize depends on the seventh       */
    /* partial derivatives of the function . since these */
    /* are inaccessible, we replace them by one in the   */
    /* estimate of deldif                                */
    /* ************************************************* */
    
    deldif   = exp(log(epsmac)/seven)*twom2;
    cgres[j] = cgres[j]+1;
    for (i = 1 ; i <= n ; i++) {
        zz    = xl[i];
        xincr = deldif*(fabs(xl[i])+one);
        xl[i] = zz-xincr;
        
        eh(j,xl,&fxl[1]);
        
        xl[i] = zz+xincr;
        
        eh(j,xl,&fxl[2]);
        
        xincr = xincr+xincr;
        d1    = xincr;
        xl[i] = zz-xincr;
        
        eh(j,xl,&fxl[3]);
        
        xl[i] = zz+xincr;
        
        eh(j,xl,&fxl[4]);
        
        xincr = xincr+xincr;
        d2    = xincr;
        xl[i] = zz-xincr;
        
        eh(j,xl,&fxl[5]);
        
        xl[i] = zz+xincr;
        
        eh(j,xl,&fxl[6]);
        
        xl[i]     = zz;
        d3        = xincr+xincr;
        sd1       = (fxl[2]-fxl[1])/d1;
        sd2       = (fxl[4]-fxl[3])/d2;
        sd3       = (fxl[6]-fxl[5])/d3;
        sd3       = sd2-sd3;
        sd2       = sd1-sd2;
        sd3       = sd2-sd3;
        gradxl[i] = sd1+p4*sd2+sd3/c45;
    }
    return;
}

/* **************************************************************************** */
/*              compute the i-th inequality constaint, bounds included          */
/* **************************************************************************** */
void eg(INTEGER i,DOUBLE x[],DOUBLE *gxi) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X
    
    cres[i+nh] = cres[i+nh]+1;
    
    /* s.t. ineq(1): */
    
    *gxi = 1.47194*x[1] - 0.312665*x[3]- 1.40314;
    
    return;
}

/* **************************************************************************** */
/*              compute the gradient of the i-th inequality constraint          */
/*          not necessary for bounds, but constant gradients must be set        */
/*                      here e.g. using dcopy from a data-field                 */
/* **************************************************************************** */
void egradg(INTEGER j,DOUBLE xl[],DOUBLE gradxl[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X
    #include "o8cons.h"
    
    void eg(INTEGER i,DOUBLE x[],DOUBLE *gxi);
    
    static DOUBLE   fxl[7];
    static DOUBLE   sd1,sd2,sd3,xincr,zz,d1,d2,d3;
    static INTEGER  i;
    
    /* ************************************************* */
    /* high precision numerical differentiation          */
    /* by sixth order extrapolation                      */
    /* uses Richardson-extrapolation of three values     */
    /* of the symmetric difference quotient with         */
    /* rather large discretization stepsize              */
    /* the optimal stepsize depends on the seventh       */
    /* partial derivatives of the function . since these */
    /* are inaccessible, we replace them by one in the   */
    /* estimate of deldif                                */
    /* ************************************************* */
    
    deldif      = exp(log(epsmac)/seven)*twom2;
    cgres[j+nh] = cgres[j+nh]+1;
    for ( i = 1 ; i <= n ; i++) {
        zz    = xl[i];
        xincr = deldif*(fabs(xl[i])+one);
        xl[i] = zz-xincr;
        
        eg(j,xl,&fxl[1]);
        
        xl[i] = zz+xincr;
        
        eg(j,xl,&fxl[2]);
        
        xincr = xincr+xincr;
        d1    = xincr;
        xl[i] = zz-xincr;
        
        eg(j,xl,&fxl[3]);
        
        xl[i] = zz+xincr;
        
        eg(j,xl,&fxl[4]);
        
        xincr = xincr+xincr;
        d2    = xincr;
        xl[i] = zz-xincr;
        
        eg(j,xl,&fxl[5]);
        
        xl[i] = zz+xincr;
        
        eg(j,xl,&fxl[6]);
        
        xl[i]     = zz;
        d3        = xincr+xincr;
        sd1       = (fxl[2]-fxl[1])/d1;
        sd2       = (fxl[4]-fxl[3])/d2;
        sd3       = (fxl[6]-fxl[5])/d3;
        sd3       = sd2-sd3;
        sd2       = sd1-sd2;
        sd3       = sd2-sd3;
        gradxl[i] = sd1+p4*sd2+sd3/c45;
    }
    return;
}

/* **************************************************************************** */
/*                        user functions (if bloc == TRUE)                      */
/* **************************************************************************** */
void eval_extern(INTEGER mode) {
    #define  X extern
    #include "o8comm.h"
    #include "o8fint.h"
    #undef   X
    #include "o8cons.h"

    return;
}