chemequi.c

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

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

main() {
    void donlp2(void);
    
    donlp2();
    
    exit(0);
}
/* Himmelblau 6 with corrections by Murtagh and Sargent */

/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
     
    static INTEGER i,j;
    
    for (i = 1 ; i <= 45 ; i++) {
        xst[i] = 0.1e0;
        x[i]   = xst[i];
    }
    strcpy(name,"himmelblau6");
    
    n    = 45;
    nh   = 16;
    ng   = 45;
    del0 = 1.e-5;
    tau0 = 1.e4;
    tau  = 0.1;
    for (j = 0 ; j <= 16 ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] = 0;
        gunit[3][j] = 0;
    }
    for (j = 17 ; j <= 61 ; j++) {
        gunit[1][j] = 1;
        gunit[3][j] = 1;
        gunit[2][j] = j-16;
    }
    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    static INTEGER i;
    
    for (i = 1 ; i <= 16 ; 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 INTEGER i,n1,n2,k,j;
    static DOUBLE  s1,s2;
    static DOUBLE  xa[NX+1];
    
    static DOUBLE  nk[] = { 0,  /* not used : index 0 */
                            4,17,35,38,41,43,45};
    
    static DOUBLE  c[]  = { 0,  /* not used : index 0 */
           0.e0    , -7.69e0   ,-11.52e0  ,-36.6e0    ,-10.94e0  ,  0.e0    ,
           0.e0    ,  0.e0     ,  0.e0    ,  0.e0     ,  0.e0    ,  2.5966e0,
         -39.39e0  ,-21.35e0   ,-32.84e0  ,  6.26e0   ,  0.e0    , 10.45e0  ,
           0.e0    , -0.5e0    ,  0.e0    ,  0.e0     ,  0.e0    ,  2.2435e0,
           0.e0    ,-39.39e0   ,-21.49e0  ,-32.84e0   ,  6.12e0  ,  0.e0    ,
           0.e0    , -1.9028e0 , -2.8889e0, -3.3622e0 , -7.4854e0,-15.639e0 ,
           0.e0    , 21.81e0   ,-16.79e0  ,  0.e0     , 18.9779e0,  0.e0    ,
           11.959e0,  0.e0     , 12.899e0};
           
    icf = icf+1;
    *fx = 0.e0;
    for (i = 1 ; i <= n ; i++) {
        xa[i] = sqrt(pow(x[i],2)+1.e-20);
    }
    for (i = 1 ; i <= 7 ; i++) {
        n1 = 1;
        n2 = nk[i];
        if ( i > 1 ) n1 = nk[i-1]+1;
        s1 = 0.e0;
        for (k = n1 ; k <= n2 ; k++) {
            s1 = s1+xa[k];
        }
        s2 = 0.e0;
        for (j = n1 ; j <= n2 ; j++) {
            s2 = s2+xa[j]*(c[j]+log(xa[j]/s1));
        }
        *fx = *fx+s2;
    }
    return;
}

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

    static INTEGER i,n1,n2,k,l,j;
    static DOUBLE  xa[NX+1],s1;
    
    static DOUBLE  nk[] = { 0,  /* not used : index 0 */
                            4,17,35,38,41,43,45};
    
    static DOUBLE  c[]  = { 0,  /* not used : index 0 */
           0.e0    , -7.69e0   ,-11.52e0  ,-36.6e0    ,-10.94e0  ,  0.e0    ,
           0.e0    ,  0.e0     ,  0.e0    ,  0.e0     ,  0.e0    ,  2.5966e0,
         -39.39e0  ,-21.35e0   ,-32.84e0  ,  6.26e0   ,  0.e0    , 10.45e0  ,
           0.e0    , -0.5e0    ,  0.e0    ,  0.e0     ,  0.e0    ,  2.2435e0,
           0.e0    ,-39.39e0   ,-21.49e0  ,-32.84e0   ,  6.12e0  ,  0.e0    ,
           0.e0    , -1.9028e0 , -2.8889e0, -3.3622e0 , -7.4854e0,-15.639e0 ,
           0.e0    , 21.81e0   ,-16.79e0  ,  0.e0     , 18.9779e0,  0.e0    ,
           11.959e0,  0.e0     , 12.899e0};

    icgf = icgf+1;
    for (i = 1 ; i <= n ; i++) {
        xa[i] = sqrt(pow(x[i],2)+1.e-20);
    }
    for (k = 1 ; k <= 7 ; k++) {
        n1 = 1;
        n2 = nk[k];
        if( k > 1 ) n1 = nk[k-1]+1;
        s1 = 0.e0;
        for (l = n1 ; l <= n2 ; l++) {
            s1 = s1+xa[l];
        }
        for (j = n1 ; j <= n2 ; j++) {
            gradf[j] =  ( c[j]+log(xa[j]/s1))*x[j]/xa[j];
        }
    }
    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;
    switch (i) {
    case 1:
        *hxi = x[ 1]+x[ 5]+x[18]+x[32]+2.e0*x[33]+3.e0*x[34]+4.e0*x[35]
                -0.6529581e0;
        break;
    case 2:
        *hxi = x[ 2]+x[ 6]+x[14]+x[15]+x[16]+x[19]+x[27]+x[28]+x[29]+x[43]
                 +x[45]-  0.281941e0;
        break;
    case 3:
        *hxi = x[ 3]+x[ 7]+x[20]-3.705233e0;
        break;
    case 4:
        *hxi = x[ 4]+x[ 8]+x[13]+x[15]-x[16]+x[21]+x[26]+x[28]-x[29]+x[36]
                 -x[38]+x[39]-x[41]-x[43]-x[45] - 47.00022e0;
        break;
    case 5:
        *hxi = x[ 4]+x[ 9]+x[13]+x[14]+x[15]+x[16]+x[22]+x[26]+x[27]+x[28]
                 +x[29]-47.02972e0;
        break;
    case 6:
        *hxi = x[10]+x[23]-0.08005e0;
        break;
    case 7:
        *hxi = x[11]+x[24]-0.08813e0;
        break;
    case 8:
        *hxi = x[12]+x[25]-0.04829e0;
        break;
    case 9:
        *hxi = x[17]-0.0155e0;
        break;
    case 10:
        *hxi = x[30]-0.0211275e0;
        break;
    case 11:
        *hxi = x[31]+x[32]+x[33]+x[34]+x[35]- 0.0022725e0;
        break;
    case 12:
        *hxi = x[ 8]-x[ 9]-x[10]+x[11]+x[12]-x[14]- 2.e0*x[16]-x[17];
        break;
    case 13:
        *hxi = +x[36]+x[37]+x[38]
                 -4.e0*x[31]-3.e0*x[32]-2.e0*x[33] -x[34];
        break;
    case 14:
        *hxi = -x[32]-2.e0*x[33]-3.e0*x[34]-4.e0*x[35]
                + x[39]+x[40]+x[41];
        break;
    case 15:
        *hxi = -4.e0*x[38]+x[42]+x[43];
        break;
    case 16:
        *hxi = -4.e0*x[41]+x[44]+x[45];
        break;
    }
    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

    static INTEGER j;

    cgres[i] = cgres[i]+1;

    for (j = 1 ; j <= 45 ; j++) {
        gradhi[j] = 0.e0;
    }
    switch (i) {
    case 1:
        gradhi[ 1] =  1.e0;
        gradhi[ 5] =  1.e0;
        gradhi[18] =  1.e0;
        gradhi[32] =  1.e0;
        gradhi[33] =  2.e0;
        gradhi[34] =  3.e0;
        gradhi[35] =  4.e0;
        
        break;
        
    case 2:
        gradhi[ 2] =  1.e0;
        gradhi[ 6] =  1.e0;
        gradhi[14] =  1.e0;
        gradhi[15] =  1.e0;
        gradhi[16] =  1.e0;
        gradhi[19] =  1.e0;
        gradhi[27] =  1.e0;
        gradhi[28] =  1.e0;
        gradhi[29] =  1.e0;
        gradhi[43] =  1.e0;
        gradhi[45] =  1.e0;
        
        break;
        
    case 3:
        gradhi[ 3] =  1.e0;
        gradhi[ 7] =  1.e0;
        gradhi[20] =  1.e0;
        
        break;
        
    case 4:
        gradhi[ 4] =  1.e0;
        gradhi[ 8] =  1.e0;
        gradhi[13] =  1.e0;
        gradhi[15] =  1.e0;
        gradhi[16] = -1.e0;
        gradhi[21] =  1.e0;
        gradhi[26] =  1.e0;
        gradhi[28] =  1.e0;
        gradhi[29] = -1.e0;
        gradhi[36] =  1.e0;
        gradhi[38] = -1.e0;
        gradhi[39] =  1.e0;
        gradhi[41] = -1.e0;
        gradhi[43] = -1.e0;
        gradhi[45] = -1.e0;
        
        break;
        
    case 5:
        gradhi[ 4] =  1.e0;
        gradhi[ 9] =  1.e0;
        gradhi[13] =  1.e0;
        gradhi[14] =  1.e0;
        gradhi[15] =  1.e0;
        gradhi[16] =  1.e0;
        gradhi[22] =  1.e0;
        gradhi[26] =  1.e0;
        gradhi[27] =  1.e0;
        gradhi[28] =  1.e0;
        gradhi[29] =  1.e0;
        
        break;
        
    case 6:
        gradhi[10] =  1.e0;
        gradhi[23] =  1.e0;
        
        break;
        
    case 7:
        gradhi[11] =  1.e0;
        gradhi[24] =  1.e0;
        
        break;
        
    case 8:
        gradhi[12] =  1.e0;
        gradhi[25] =  1.e0;
        
        break;
        
    case 9:
        gradhi[17] =  1.e0;
        
        break;
        
    case 10:
        gradhi[30] =  1.e0;
        
        break;
        
    case 11:
        gradhi[31] =  1.e0;
        gradhi[32] =  1.e0;
        gradhi[33] =  1.e0;
        gradhi[34] =  1.e0;
        gradhi[35] =  1.e0;
        
        break;
        
    case 12:
        gradhi[ 8] =  1.e0;
        gradhi[ 9] = -1.e0;
        gradhi[10] = -1.e0;
        gradhi[11] =  1.e0;
        gradhi[12] =  1.e0;
        gradhi[14] = -1.e0;
        gradhi[16] = -2.e0;
        gradhi[17] = -1.e0;
        
        break;
        
    case 13:
        gradhi[31] = -4.e0;
        gradhi[32] = -3.e0;
        gradhi[33] = -2.e0;
        gradhi[34] = -1.e0;
        gradhi[36] =  1.e0;
        gradhi[37] =  1.e0;
        gradhi[38] =  1.e0;
        
        break;
        
    case 14:
        gradhi[32] = -1.e0;
        gradhi[33] = -2.e0;
        gradhi[34] = -3.e0;
        gradhi[35] = -4.e0;
        gradhi[39] =  1.e0;
        gradhi[40] =  1.e0;
        gradhi[41] =  1.e0;
        
        break;
        
    case 15:
        gradhi[38] = -4.e0;
        gradhi[42] =  1.e0;
        gradhi[43] =  1.e0;
        
        break;
        
    case 16:
        gradhi[41] = -4.e0;
        gradhi[44] =  1.e0;
        gradhi[45] =  1.e0;
        
        break;
    }
    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
    
    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;
    
    *gxi = x[i]-1.e-6;
    
    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;
    for (j = 1 ; j <= 45 ; j++) {
        gradgi[j] = 0.e0;
    }
    gradgi[i] =  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;
}