him13.c

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

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


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

/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {

    /*  this is example Himmelblau 13, box parameter identification problem */

    #define  X extern
    #include "o8comm.h"
    #undef   X
    void startx(DOUBLE x[]);    
    static INTEGER i,j;

    strcpy(name,"him13boxproblem");

    n  = 12;
    nh = 7;
    ng = 16;
    
    del0 = 0.1e0;
    tau0 = 1.0e0;
    tau  =  .1e0;
    
    x[1] = 2.52e0;
    x[2] = 2.e0;
    x[3] = 37.5e0;
    x[4] = 9.25e0;
    x[5] = 6.8e0;
    
    for (i = 6 ; i <= n ; i++) {
        x[i] = 0.0e0;
    }
    startx(x);
    for (j = 0 ; j <= 7 ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] = 0;
        gunit[3][j] = 0;
    }
    for (i = 8 ; i <= 15 ; i++) {
        gunit[1][j] = 1;
        gunit[2][i] = i - 7;
        gunit[3][i] = 1;
    }
    for (i = 16 ; i <= 23 ; i++) {
        gunit[1][j] = 1;
        gunit[2][i] = i - 15;
        gunit[3][i] = -1;
    }
    return;
}

void startx(DOUBLE x[]) {
    static INTEGER  i;
    static DOUBLE   xk[31] = {0.,/* not used : index 0 */
            -145.421402e3 ,   2.9311506e3,-40.427932e0  ,  5.106192e3,
              15.711360e3 ,-161.622577e3 ,  4.17615328e3,  2.8260078e0,
               9.200476e3 ,  13.160295e3 ,-21.6869194e3 ,123.56928e0,
             -21.1188894e0, 706.834000e0 ,  2.8985730e3 , 28.298388e3,
              60.8109600e0,  31.242116e0 ,329.5740000e0 , -2.882082e3,
              74.0953845e3, -30.6262544e1, 16.2436490e0 , -3.094252e3,
             - 5.5662628e3, -26.237000e3 , 99.e0        , -0.42e0,  
               1.3e3      ,   2.1e3};
            
    x[6] = xk[1];
    for ( i = 2 ; i <= 5 ; i++) {
        x[6] = x[6]+xk[i]*x[i];
    }
    x[6] = x[6]*x[1];
    x[9] = xk[6];
    for ( i = 2 ; i <= 5 ; i++) {
        x[9] = x[9]+xk[i+5]*x[i];
    }
    x[10] = xk[11];
    for ( i = 2 ; i <= 5 ; i++) {
        x[10] = x[10]+xk[i+10]*x[i];
    }
    x[11] = xk[16];
    for ( i = 2 ; i <= 5 ; i++) {
        x[11] = x[11]+xk[i+15]*x[i];
    }
    x[12] = xk[21];
    for ( i = 2 ; i <= 5 ; i++) {
        x[12] = x[12]+xk[i+20]*x[i];
    }
    x[7] = (x[9]+x[10]+x[11])*x[1];
    x[8] = xk[26];
    for ( i = 2 ; i <= 5 ; i++) {
        x[8] = x[8]+xk[i+25]*x[i];
    }
    x[8] = x[8]*x[1]+x[6]+x[7];
    
    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X

    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 xk[6] ={0.,/* not used : index 0 */
    925.548252e3,-61.9688432e3,23.3088196e0,-27.097648e3,-50.843766e3};
    
    icf = icf+1;
    
    
    *fx = -(50.e0*x[9]+9.583e0*x[10]+20.e0*x[11]+15.e0*x[12]-
             852.96e3-38.1e3*(x[2]+0.01e0*x[3]) +xk[1]+xk[2]*x[2]+
             xk[3]*x[3]+xk[4]*x[4]+xk[5]*x[5])*x[1]-15.e0*x[6];
    
    return;
}

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

    static DOUBLE xk[6] = {0.,/* not used : index 0 */
    925.548252e3,-61.9688432e3,23.3088196e0,-27.097648e3,-50.843766e3};
    
    icgf = icgf+1;
    
    gradf[1]  = -(50.e0*x[9]+9.583e0*x[10]+20.e0*x[11]+15.e0*x[12]-
                852.96e3-38.1e3*(x[2]+0.01e0*x[3]) +xk[1]+xk[2]*x[2]+
                xk[3]*x[3]+xk[4]*x[4]+xk[5]*x[5]);
    gradf[2]  = 38.1e3*x[1]-xk[2]*x[1];
    gradf[3]  = -xk[3]*x[1] +.381e3*x[1];
    gradf[4]  = -xk[4]*x[1];
    gradf[5]  = -xk[5]*x[1];
    gradf[6]  = -15.e0;
    gradf[7]  =  0.0e0;
    gradf[8]  =  0.0e0;
    gradf[9]  = -x[1]*50.e0;
    gradf[10] = -x[1]*9.583e0;
    gradf[11] = -x[1]*20.e0;
    gradf[12] = -x[1]*15.e0;
    
    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

    
    static INTEGER k,ii,iik;
    static DOUBLE xk[31] = {0.,/* not used : index 0 */
            -145.421402e3 ,   2.9311506e3,-40.427932e0  ,  5.106192e3,
              15.711360e3 ,-161.622577e3 ,  4.17615328e3,  2.8260078e0,
               9.200476e3 ,  13.160295e3 ,-21.6869194e3 ,123.56928e0,
             -21.1188894e0, 706.834000e0 ,  2.8985730e3 , 28.298388e3,
              60.8109600e0,  31.242116e0 ,329.5740000e0 , -2.882082e3,
              74.0953845e3, -30.6262544e1, 16.2436490e0 , -3.094252e3,
             - 5.5662628e3, -26.237000e3 , 99.e0        , -0.42e0,
               1.3e3      ,   2.1e3};
               
    cres[i] = cres[i]+1;
    
    
    if (i != 7) {
        ii   = (i-1)*5;
        *hxi = xk[ii+1];
        for ( k = 2 ; k <= 5 ; k++) {
            iik  = ii+k;
            *hxi = *hxi+xk[iik]*x[k];
        }
        switch (i) {
        case 1:
            *hxi = x[6]-*hxi*x[1];
            break;
        case 2:
        case 3:
        case 4:
        case 5:
            *hxi = x[i+7]-*hxi;
            break;
        case 6:
            *hxi = x[8]-(*hxi*x[1]+x[6]+x[7]);
            break;
        }
    } else {
        *hxi = x[7]-(x[9]+x[10]+x[11])*x[1];
    }
    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  k,ii,iik;
    static DOUBLE   xk[31] = {0.,/* not used : index 0 */
            -145.421402e3 ,   2.9311506e3,-40.427932e0  ,  5.106192e3,
              15.711360e3 ,-161.622577e3 ,  4.17615328e3,  2.8260078e0,
               9.200476e3 ,  13.160295e3 ,-21.6869194e3 ,123.56928e0,
             -21.1188894e0, 706.834000e0 ,  2.8985730e3 , 28.298388e3,
              60.8109600e0,  31.242116e0 ,329.5740000e0 , -2.882082e3,
              74.0953845e3, -30.6262544e1, 16.2436490e0 , -3.094252e3,
              -5.5662628e3, -26.237000e3 , 99.e0        , -0.42e0,
               1.3e3      ,   2.1e3};
               
    cgres[i] = cgres[i]+1;
    for ( k = 1 ; k <= 12 ; k++) {
        gradhi[k] = 0.e0;
    }
    switch (i) {
    case 1:
    case 6:
        ii = (i-1)*5;
        gradhi[1] = -(xk[ii+1]+xk[ii+2]*x[2]+xk[ii+3]*x[3]+xk[ii+4]*x[4]
                    +xk[ii+5]*x[5]);
        for ( k = 2 ; k <= 5 ; k++) {
            iik       = ii+k;
            gradhi[k] = -xk[iik]*x[1];
        }
        if (i == 6) {
            gradhi[6] = -1.e0;
            gradhi[7] = -1.e0;
            gradhi[8] =  1.e0;
        } else {
            gradhi[6] = 1.e0;
        }
        break;
    case 2:
    case 3:
    case 4:
    case 5:
        ii = (i-1)*5;
        for ( k = 2 ; k <= 5 ; k++) {
            iik = ii+k;
            gradhi[k] = -xk[iik];
        }
        gradhi[i+7] = 1.e0;
        break;
    case 7:
        gradhi[ 1] = -(x[9]+x[10] +x[11]);
        gradhi[ 7] =  1.e0;
        gradhi[ 9] = -x[1];
        gradhi[10] = -x[1];
        gradhi[11] = -x[1];
        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 DOUBLE xk[17] = {0.,/* not used : index 0 */
            0.e0,1.2e0,20.0e0,9.e0 , 6.5e0,  0.e0,  0.e0,  0.e0,
            5.e0,2.4e0, 6.e1 ,9.3e0, 7.0e0,294.e3,294.e3,277.2e3};
            
    if(i <= 8) {
        *gxi = x[i]-xk[i];
    } else {
        *gxi = xk[i]-x[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;
    
    for ( j = 1 ; j <= 12 ; j++) {
        gradgi[j] = 0.e0;
    }
    if(i <= 8) {
        gradgi[i]   = 1.0e0;
    } else {
        gradgi[i-8] = -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;
}