alkylati.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,j;
    static DOUBLE   xst0[11] = {0.,/* not used : index 0 */
                                1745.e0  ,12.e3  ,11.e1,3048.e0 ,1974.e0,
                                  89.2e0 ,92.8e0 , 8.e0,   3.6e0, 145.e0 };
                                      
    /* 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,"alkylation");
     
    /* x is initial guess and also holds the current solution */
    /* problem dimension n = dim(x), nh = dim(h), ng = dim(g) */
    
    n      = 10;
    nh     = 3;
    ng     = 28;
    analyt = TRUE;
    epsdif = 1.e-16;
    
    /* epsfcn   = 1.e-16; */
    /* taubnd   = 5.e-6; */
    /*  bloc    = TRUE; */
    /* difftype = 3; */
    
    nreset = n;
    
    /* del0 and tau0: see below */
    
    del0 = 0.2e0;
    tau0 = 1.e0;
    tau  = 0.1e0;
    for (i = 1 ; i <= n ; i++) {
        x[i] = xst0[i];
    }
    /* gunit-array, see donlp2doc.txt */
    
    for (j = 0 ; j <= 11 ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] = 0;
        gunit[3][j] = 0;
    }
    for (j = 12 ; j <= 31 ; j++) {
        gunit[1][j] = 1;
        if ( j <= 21 ) {
            gunit[2][j] = j-11;
            gunit[3][j] = 1;
        } else {
            gunit[2][j] = j-21;
            gunit[3][j] = -1;
        }
    }
    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    te2=TRUE;
    te3=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

    icf = icf+1;
    *fx = 5.04e0*x[1]+.035e0*x[2]+10.e0*x[3]+3.36e0*x[5]-.063e0*x[4]*x[7];
    
    return;
}

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

    static INTEGER  j;
    static DOUBLE   a[11] = {0.,/* not used : index 0 */
                             5.04e0,0.035e0,10.e0,0.e0,3.36e0,
                             0.e0  ,0.e0   , 0.e0,0.e0,0.e0};

    icgf = icgf+1;
    for (j = 1 ; j <= 10 ; j++) {
        gradf[j] = a[j];
    }
    gradf[4] = -0.063e0*x[7];
    gradf[7] = -0.063e0*x[4];
    
    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 = 1.22e0*x[4]-x[1]-x[5];
        break;
    case 2:
        *hxi = 9.8e4*x[3]/(x[4]*x[9]+1.e3*x[3])-x[6];
        break;
    case 3:
        *hxi = (x[2]+x[5])/x[1]-x[8];
        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;
    static DOUBLE   t,t1;
    
    cgres[i] = cgres[i]+1;
    for (j = 1 ; j <= 10 ; j++) {
        gradhi[j] = 0.e0;
    }
    switch (i) {
    case 1:
        gradhi[1] = -1.e0;
        gradhi[4] = 1.22e0;
        gradhi[5] = -1.e0;
        break;
    case 2:
        t         = 9.8e4/(x[4]*x[9]+1.e3*x[3]);
        t1        = t/(x[4]*x[9]+1.e3*x[3])*x[3];
        gradhi[3] = t-1.e3*t1;
        gradhi[4] = -x[9]*t1;
        gradhi[9] = -x[4]*t1;
        gradhi[6] = -1.e0;
        break;
    case 3:
        gradhi[1] = -(x[2]+x[5])/pow(x[1],2);
        gradhi[2] = 1.e0/x[1];
        gradhi[5] = gradhi[2];
        gradhi[8] = -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
    
    static INTEGER  k;
    static DOUBLE   og[11] = {0.,/* not used : index 0 */
                              2.e3 ,16.e3, 1.2e2,5.e3,  2.e3,
                              93.e0,95.e0,12.e0, 4.e0,162.e0 };

    static DOUBLE   ug[11] = {0.,/* not used : index 0 */
                              1.e-5, 1.e-5,1.e-5,1.e-5,  1.e-5,
                             85.e0, 90.e0 ,3.e0 ,1.2e0,145.e0 };
    
    static DOUBLE   a = .99e0,b = .9e0,c = 2.01010101010101e-2,
                    d = 2.11111111111111e-1;
    static DOUBLE   t;
    
    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;
    
    k = (i+1)/2;

    switch (k) {
    case 1:
        t = 35.82e0-.222e0*x[10]-b*x[9];
        if(k+k == i) t = -t+x[9]*d;
        *gxi = t;
        break;
    case 2:
        t = -133.e0+3.e0*x[7]-a*x[10];
        if(k+k == i) t = -t+c*x[10];
        *gxi = t;
        break;
    case 3:
        t = 1.12e0*x[1]+.13167e0*x[1]*x[8]-.00667e0*x[1]*pow(x[8],2)-a*x[4];
        if(k+k == i) t = -t+c*x[4];
        *gxi = t;
        break;
    case 4:
        t = 57.425e0+1.098e0*x[8]-.038e0*pow(x[8],2)+.325e0*x[6]-a*x[7];
        if(k+k == i) t = -t+c*x[7];
        *gxi = t;
        break;
    default:
        if( i  > 18 ) *gxi = og[i-18]-x[i-18];
        if( i <= 18 ) *gxi = x[i-8]-ug[i-8];
    }
    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,k;
    static DOUBLE   a = .99e0,b = .9e0,c = 1.01010101010101e0,
                    d = 1.11111111111111e0;

    for ( j = 1 ; j <= NX ; j++) {
         gradgi[j] = 0.e0;
    }
    k = (i+1)/2;
    
    switch (k) {
    case 1:
        if ( k+k != i ) {
            gradgi[ 9] = -b;
            gradgi[10] = -.222e0;
        } else {
            gradgi[ 9] =  d;
            gradgi[10] =  .222e0;
        }
        break;
    case 2:
        if ( k+k != i ) {
            gradgi[ 7] =  3.e0;
            gradgi[10] = -a;
        } else {
            gradgi[ 7] = -3.e0;
            gradgi[10] =  c;
        }
        break;
    case 3:
        gradgi[1] = 1.12e0+.13167e0*x[8]-.00667e0*pow(x[8],2);
        gradgi[4] = -a;
        gradgi[8] = .13167e0*x[1]-.01334e0*x[1]*x[8];
        if( k+k == i ) {
            gradgi[1] = -gradgi[1];
            gradgi[8] = -gradgi[8];
            gradgi[4] = c;
        }
        break;
    case 4:
        gradgi[6] = .325e0;
        gradgi[7] = -a;
        gradgi[8] = 1.098e0-.076e0*x[8];
        if(k+k == i) {
            gradgi[6] = -.325e0;
            gradgi[7] = c;
            gradgi[8] = -gradgi[8];
        }
        break;
    default:
        if( i  > 18 ) gradgi[i-18] = -1.e0;
        if( i <= 18 ) gradgi[i-8]  = 1.e0;
        break;
    }
    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;
}