robota.c

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

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

main() {
    void donlp2(void);
    
    donlp2();
    
    exit(0);
}
/* coded by: Hans D. Mittelmann <[e-mail]> (*) */

#define XDIM    21
#define BDM     81

static DOUBLE  c11 = 2.e0, c12 =  8.e0, c13 =  250.e0;
static DOUBLE  c21 = 3.e0, c22 = 18.e0, c23 =  650.e0;
static DOUBLE  c31 = 4.e0, c32 = 50.e0, c33 = 1000.e0;

static DOUBLE  vi[BDM+1][XDIM+1],vi1[BDM+1][XDIM+1],vi2[BDM+1][XDIM+1],
               v11[BDM+1],v12[BDM+1],v13[BDM+1],
               v21[BDM+1],v22[BDM+1],v23[BDM+1],
               v31[BDM+1],v32[BDM+1],v33[BDM+1];
    
static DOUBLE  bb[XDIM+1];

/* (*) translation of -, from f77 to C by S. Schoeffert, ASB, Bourges,France    */
/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    DOUBLE th11 (DOUBLE t);
    DOUBLE th12 (DOUBLE t);
    DOUBLE th13 (DOUBLE t);
    DOUBLE th21 (DOUBLE t);
    DOUBLE th22 (DOUBLE t);
    DOUBLE th23 (DOUBLE t);
    DOUBLE th31 (DOUBLE t);
    DOUBLE th32 (DOUBLE t);
    DOUBLE th33 (DOUBLE t);
    DOUBLE phil (INTEGER i, DOUBLE t);
    DOUBLE phil1(INTEGER i, DOUBLE t);
    DOUBLE phil2(INTEGER i, DOUBLE t);

    static INTEGER i,j,m;
    static DOUBLE  xst0[XDIM+1],h,tp;

    /*              other data declaration for cases a and b                    */
    /*                                                                          */
    /* a    for (i = 1 ; i <= XDIM ; i++) xst0[i] = 1.491400623321533e0;        */
    /* b    for (i = 1 ; i <= XDIM ; i++) xst0[i] = 4.377146244049071e0;        */
    /* c    for (i = 1 ; i <= XDIM ; i++) xst0[i] = 1.920426464080810e0;        */
    /* a    static DOUBLE  c11 = 2.e0, c12 =  8.e0, c13 =  250.e0;              */
    /* a    static DOUBLE  c21 = 3.e0, c22 = 18.e0, c23 =  650.e0;              */
    /* a    static DOUBLE  c31 = 4.e0, c32 = 50.e0, c33 = 1000.e0;              */
    /* b    static DOUBLE  c11 = 1.e0, c12 =  3.e0, c13 =  100.e0;              */
    /* b    static DOUBLE  c21 = 1.e0, c22 =  3.e0, c23 =  100.e0;              */
    /* b    static DOUBLE  c31 = 1.e0, c32 =  3.e0, c33 =  100.e0;              */
    /* c    static DOUBLE  c11 = 2.e0, c12 =  8.e0, c13 =   25.e0;              */
    /* c    static DOUBLE  c21 = 3.e0, c22 = 18.e0, c23 =   65.e0;              */
    /* c    static DOUBLE  c31 = 4.e0, c32 = 50.e0, c33 =  100.e0;              */
    /* a    strcpy(name,"robota");                                              */
    /* b    strcpy(name,"robotb");                                              */

    n = XDIM;
    m = BDM;
    for (i = 1 ; i <= XDIM ; i++) {
        xst0[i] = 1.491400623321533e0;
    }
    strcpy(name,"robota");

    nh     = 0;
    ng     = m*18;
    h      = 1.e0/(n-1);
    silent = FALSE;
    analyt = TRUE;
    epsdif = 0.e0;
    del0   = 1.e0;
    tau0   = 1.e4;
    for (j = 1 ; j <= m ; j++) {
        tp     = (DOUBLE)(j-1)/(DOUBLE)(m-1);
        v11[j] = th11(tp);
        v12[j] = th12(tp);
        v13[j] = th13(tp);
        v21[j] = th21(tp);
        v22[j] = th22(tp);
        v23[j] = th23(tp);
        v31[j] = th31(tp);
        v32[j] = th32(tp);
        v33[j] = th33(tp);
        for (i = 1 ; i <= n ; i++) {
            vi [j][i] = phil(i,tp);
            vi1[j][i] = phil1(i,tp);
            vi2[j][i] = phil2(i,tp);
        }
    }
    for (i = 1 ; i <= n ; i++) {
        x[i]  = xst0[i];
        bb[i] = h;
    }
    bb[1]   = .5e0*h;
    bb[n]   = bb[1];
    bb[2]   = 23.e0*h/24.e0;
    bb[n-1] = bb[2];
    for (j = 0 ; j <= nh+ng ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] = 0;
        gunit[3][j] = 0;
    }
    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    /* change termination criterion */
    
    epsx = 1.e-7;
    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;
    
    icf = icf+1;
    *fx = 0.e0;
    for (i = 1 ; i <= n ; i++) {
        *fx = *fx+bb[i]*x[i];
    }
    return;
}

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

    static INTEGER i;
    
    icgf = icgf+1;
    for (i = 1 ; i <= n ; i++) {
        gradf[i] = bb[i];
    }
    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;

    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

    if ( gunit[1][i] != 1 ) cgres[i] = cgres[i]+1;
    
    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 ip,k,j;
    static DOUBLE  sum,sum1,sum2;

    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;
    
    ip   = (i-1)/18+1;
    k    = i-18*(ip-1);
    sum  = 0.e0;
    sum1 = 0.e0;
    sum2 = 0.e0;
    for (j = 1 ; j <= n ; j++) {
        sum  = sum +vi [ip][j]*x[j];
        sum1 = sum1+vi1[ip][j]*x[j];
        sum2 = sum2+vi2[ip][j]*x[j];
    }
    switch (k) {
    case 1:
        *gxi = c11*sum-v11[ip];
        return;
    case 2:
        *gxi = c11*sum+v11[ip];
        return;
    case 3:
        *gxi = c12*pow(sum,3)-(v12[ip]*sum-v11[ip]*sum1);
        return;
    case 4:
        *gxi = c12*pow(sum,3)+(v12[ip]*sum-v11[ip]*sum1);
        return;
    case 5:
        *gxi = c13*pow(sum,5)-(v13[ip]*pow(sum,2)-3.e0*v12[ip]*sum*sum1
               +3.e0*v11[ip]*pow(sum1,2)+v11[ip]*sum*sum2);
        return;
    case 6:
        *gxi = c13*pow(sum,5)+(v13[ip]*pow(sum,2)-3.e0*v12[ip]*sum*sum1
               +3.e0*v11[ip]*pow(sum1,2)+v11[ip]*sum*sum2);
        return;
    case 7:
        *gxi = c21*sum-v21[ip];
        return;
    case 8:
        *gxi = c21*sum+v21[ip];
        return;
    case 9:
        *gxi = c22*pow(sum,3)-(v22[ip]*sum-v21[ip]*sum1);
        return;
    case 10:
        *gxi = c22*pow(sum,3)+(v22[ip]*sum-v21[ip]*sum1);
        return;
    case 11:
        *gxi = c23*pow(sum,5)-(v23[ip]*pow(sum,2)-3.e0*v22[ip]*sum*sum1
               +3.e0*v21[ip]*pow(sum1,2)+v21[ip]*sum*sum2);
        return;
    case 12:
        *gxi = c23*pow(sum,5)+(v23[ip]*pow(sum,2)-3.e0*v22[ip]*sum*sum1
               +3.e0*v21[ip]*pow(sum1,2)+v21[ip]*sum*sum2);
        return;
    case 13:
        *gxi = c31*sum-v31[ip];
        return;
    case 14:
        *gxi = c31*sum+v31[ip];
        return;
    case 15:
        *gxi = c32*pow(sum,3)-(v32[ip]*sum-v31[ip]*sum1);
        return;
    case 16:
        *gxi = c32*pow(sum,3)+(v32[ip]*sum-v31[ip]*sum1);
        return;
    case 17:
        *gxi = c33*pow(sum,5)-(v33[ip]*pow(sum,2)-3.e0*v32[ip]*sum*sum1
               +3.e0*v31[ip]*pow(sum1,2)+v31[ip]*sum*sum2);
        return;
    case 18:
        *gxi = c33*pow(sum,5)+(v33[ip]*pow(sum,2)-3.e0*v32[ip]*sum*sum1
               +3.e0*v31[ip]*pow(sum1,2)+v31[ip]*sum*sum2);
        return;
    }
}

/* **************************************************************************** */