weaponas.c

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

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

main() {
    void donlp2(void);
    
    donlp2();
    
    exit(0);
}
static DOUBLE  wd[101];

static DOUBLE  wa[][21] = {
{0.,    0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,
            0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    ,0.    },
        
{0.,    1.e0  , .95e0,1.e0  ,1.e0  ,1.e0  , .85e0, .9e0 , .85e0, .8e0 ,1.e0  ,  
            1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  , .95e0,1.e0  ,1.e0  }, 
{0.,     .84e0, .83e0, .85e0, .84e0, .85e0, .81e0, .81e0, .82e0, .8e0 , .86e0,  
            1.e0  , .98e0,1.e0  , .88e0, .87e0, .88e0, .85e0, .84e0, .85e0, .85e0}, 
{0.,     .96e0, .95e0, .96e0, .96e0, .96e0, .90e0, .92e0, .91e0, .92e0, .95e0,  
             .99e0, .98e0, .99e0, .98e0, .97e0, .98e0, .95e0, .92e0, .93e0, .92e0}, 
{0.,    1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  ,1.e0  , .96e0,  
             .91e0, .92e0, .91e0, .92e0, .98e0, .93e0,1.e0  ,1.e0  ,1.e0  ,1.e0  }, 
{0.,     .92e0, .94e0, .92e0, .95e0, .95e0, .98e0, .98e0,1.e0  ,1.e0  , .90e0,
             .95e0, .96e0, .91e0, .98e0, .99e0, .99e0,1.e0  ,1.e0  ,1.e0  ,1.e0  }};
    
static DOUBLE  wc[] = {0., /* not used : index 0 */
                    2.e2,1.e2,3.e2,1.5e2,2.5e2};

static DOUBLE  wb[] = {0., /* not used : index 0 */
                    30.e0,1.e2,4.e1,5.e1,7.e1,35.e0,1.e1};

static INTEGER wj[] = {0, /* not used : index 0 */
                    1,6,10,14,15,16,20};

static DOUBLE  wu[] = {0., /* not used : index 0 */
                    6.e1,5.e1,5.e1,75.e0,4.e1,6.e1,3.5e1,3.e1,2.5e1,1.5e2,3.e1,
                    4.5e1,1.25e2,2.e2,2.e2,1.3e2,1.e2,1.e2,1.e2,1.5e2};

/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    static INTEGER i,j,nuser = 100;

    static DOUBLE xst0[] = {0., /* not used : index 0 */
         6.e0, 6.e0, 6.e0, 6.e0, 6.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0,20.e0,20.e0,20.e0,20.e0,20.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 8.e0, 8.e0, 8.e0, 8.e0, 8.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0,10.e0,10.e0,10.e0,10.e0,10.e0,
        12.e0,12.e0,12.e0,12.e0,12.e0, 7.e0, 7.e0, 7.e0, 7.e0, 7.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 1.e0,
         1.e0, 1.e0, 1.e0, 1.e0, 1.e0, 3.e0, 3.e0, 3.e0, 3.e0, 3.e0  };
        
    /* name is ident of the example/user and can be set at users will           */
    /* the first static char [] must be alphabetic.  40 static char [s] maximum */

    strcpy(name,"weaponh23");

    /* x is initial guess and also holds the current solution */
    /* problem dimension n = dim(x), nh = dim(h), ng = dim(g) */

    n      = nuser;
    nh     = 0;
    ng     = 112;
    analyt = TRUE;
    epsdif = 0;
    del0   = 1.e1;
    tau0   = 1.e4;
    tau    = 0.1e0;
    for (i = 1 ; i <= n ; i++) {
        x[i] = xst0[i];
    }
    for (j = 0 ; j <= 12 ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] = 0;
        gunit[3][j] = 0;
    }
    for (j = 13 ; j <= 112 ; j++) {
        gunit[1][j] = 1;
        gunit[2][j] = j-12;
        gunit[3][j] = 1;
    }
    for (i = 1 ; i <= 5 ; i++) {
        for (j = 1 ; j <= 20 ; j++) {
            wd[(i-1)*20+j] = log(wa[i][j]);
        }
    }
    nreset = n;
    
    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    static INTEGER i;
    
    gconst[0] = FALSE;
    for (i = 1 ; i <= 112 ; i++) {
        gconst[i] = TRUE;
    }
    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

    static DOUBLE  sum,prod;
    static INTEGER i,j,k,ind;

    icf = icf+1;
    sum = 0.e0;
    for (j = 1 ; j <= 20 ; j++) {
        prod = 0.e0;
        k    = j-20;
        for (i = 1 ; i <= 5 ; i++) {
            ind  = i*20+k;
            prod = prod+x[ind]*wd[ind];
        }
        prod = exp(prod)-1.e0;
        sum  = sum+wu[j]*prod;
    }
    *fx = sum;
    
    return;
}

/* **************************************************************************** */
/*                          gradient of objective function                      */
/* **************************************************************************** */
void egradf(DOUBLE x[],DOUBLE gradf[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X

    static DOUBLE  sum;
    static INTEGER i,j,k,ind;

    icgf = icgf+1;
    for (j = 1 ; j <= 20 ; j++) {
        sum = 0.e0;
        k   = j-20;
        for (i = 1 ; i <= 5 ; i++) {
            ind = i*20+k;
            sum = sum+x[ind]*wd[ind];
        }
        for (i = 1 ; i <= 5 ; i++) {
            ind        = i*20+k;
            gradf[ind] = wu[j]*wd[ind]*exp(sum);
        }
    }
    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
    
    return;
}

/* **************************************************************************** */
/*              compute the gradient of the i-th equality constraint            */
/* **************************************************************************** */
void egradh(INTEGER i,DOUBLE x[],DOUBLE gradhi[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X

    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
    
    static DOUBLE  sum;
    static INTEGER i1,i2,j;

    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;
    
    if ( i  > 12 ) {
        *gxi = x[i-12];
    } else if ( i  > 7 ) {
        sum = 0.e0;
        i1  = i-7;
        i2  = (i1-1)*20;
        for (j = 1 ; j <= 20 ; j++) {
            sum = sum+x[i2+j];
        }
        *gxi = -sum+wc[i1];
    } else {
        sum = 0.e0;
        for (j = 1 ; j <= 5 ; j++) {
            sum = sum+x[(j-1)*20+wj[i]];
        }
        *gxi = sum-wb[i];
    }
    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 i,DOUBLE x[],DOUBLE gradgi[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X

    static INTEGER j,i1,i2;

    if ( gunit[1][i+nh] != 1 ) cgres[i+nh] = cgres[i+nh]+1;
    if ( gunit[1][i+nh] == 1 ) return;
    
    for ( j = 1 ;  j <= NX ;  j++) {
        gradgi[j] = 0.e0;
    }
    if ( i <= 7 ) {
        for (j = 1 ; j <= 5 ; j++) {
            gradgi[(j-1)*20+wj[i]] = 1.e0;
        }
    } else {
        i1 = i-7;
        i2 = (i1-1)*20;
        for (j = 1 ; j <= 20 ; j++) {
            gradgi[i2+j] = -1.e0;
        }
    }
    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;
}