paul_brown.c

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

/* ****************************************************************************
*/
/*                                 user functions
*/
/* ****************************************************************************
*/
#include "o8para.h"
#define PROP   1
#define MU1    2
#define SD1    3
#define MU2    4
#define SD2    5
main() {
    void donlp2(void);

    donlp2();

    exit(0);
}
/*  paul brown problem           */
/*          */

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

    static INTEGER i,j,nuser = 5;
    static DOUBLE  xst0[] = {0., /* not used : index 0 */
  0.0,0.1,0.1,0.1,0.1 };

/* 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,"paulbrown");

    /* 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     = 12;
    analyt = TRUE;
    epsdif = 1.e-16;
/*  old    del0   = 1.e0; */
    del0 = 1.e0;
    tau0   = 1.e1;

    for (i = 1 ; i <= n ; i++) {
        x[i] = xst0[i];
    }
    /* gunit-array, see donlpdoc.ps */

     for (i=0;i<=2;i++)
     {
       gunit[1][i] = -1;
       gunit[2][i] = 0; 
       gunit[3][i] = 0; 
     }
      
     /* this implies that in eg/egradg the two general constraints
        g1 and g2 come first */
        
     for ( i=3; i<=12; i++ )
     {
        gunit[1][i] = 1 ;
        gunit[2][i] = ( i<=7 )? i-2 : i-7 ;  /* 1,2,3,4,5,1,2,3,4,5 */
        gunit[3][i] = ( i<=7 )? 1  :  -1 ;  /* lower bounds come first */
     }


    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


    *fx = 0.05*(pow(x[1],2)+pow(x[2],2)+pow(x[3],2)+pow(x[4],2))
             +0.025*pow(x[5],2);
    
    return;
}

/* **************************************************************************** */
/*                          gradient of objective function                      */
/* **************************************************************************** */
void egradf(DOUBLE x[],DOUBLE gradf[]) {
    #define  X extern
    #include "o8fuco.h"
    #undef   X
     gradf[PROP]= 0.1  * x[PROP];
     gradf[MU1] = 0.1  * x[MU1];
     gradf[SD1] = 0.1  * x[SD1];
     gradf[MU2] = 0.1  * x[MU2];
     gradf[SD2] = 0.05 * x[SD2];


    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   og[6] = 
    {0.0,0.9999,50.0,99.0,50.0,99.0};
    static DOUBLE   ug[6] = 
    {0.0,0.0001,1.0e-6,0.0000000001,1.0e-6,0.0000000001};
   
    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;
    switch (i) {

    case 1:
      *gxi = (x[MU1] + pow(x[SD1],2)/2.0)
           - (x[MU2] + pow(x[SD2],2)/2.0);
      break;
    
    case 2:
      if ( 2.0 * x[MU1] + pow(x[SD1],2) <= 0.0 ||
           2.0 * x[MU2] + pow(x[SD2],2)<= 0.0  )
      {
        cfuerr[2] = TRUE ;
        *gxi = 0 ;
        break ;
      }
      *gxi = log(2.0 * x[MU1] + pow(x[SD1],2)) * (log(pow(x[SD1],2)) -1.0)
           - log(2.0 * x[MU2] + pow(x[SD2],2)) * (log(pow(x[SD2],2)) -1.0);
      break;
    
    default:
      if (i > 7) *gxi = og[i-7] -  x[i-7];
      else       *gxi =  x[i-2] - ug[i-2];
    }   




}

/* **************************************************************************** */
/*              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
    
    INTEGER j;   
    if ( i == 1 )
    {
      cgres[1] += 1;
      gradgi[PROP] = 0.0;   
      gradgi[MU1] = 1.0 ;   
      gradgi[SD1] = x[SD1] ;
      gradgi[MU2] = -1.0;
      gradgi[SD2] = -x[SD2] ;
    }
    if ( i ==2 )
    {
      cgres[2] += 1;
      gradgi[PROP] = 0.0;
      gradgi[MU1] = 2.0/(2.0 * x[MU1] + pow(x[SD1],2))* 
              (log(pow(x[SD1],2)) -1.0);
      gradgi[SD1] = 2.0*x[SD1]/(2.0 * x[MU1] + pow(x[SD1],2) )
                               * (log(pow(x[SD1],2)) -1.0)
                  +log(2.0 * x[MU1] + pow(x[SD1],2)) *2.0/x[SD1];
      gradgi[MU2] = 2.0/(2.0 * x[MU2] + pow(x[SD2],2))* 
              (log(pow(x[SD2],2)) -1.0);
      gradgi[SD2] = 2.0*x[SD2]/(2.0 * x[MU2] + pow(x[SD2],2) )
                               * (log(pow(x[SD2],2)) -1.0)
                  +log(2.0 * x[MU2] + pow(x[SD2],2) ) *2.0/x[SD2];
                  
    }        /* no gradients for the bounds need to be provided , these are computed from
    gunit */
    return; 
}
 




void eval_extern(INTEGER mode) {
    #define  X extern
    #include "o8comm.h"
    #include "o8fint.h"
    #undef   X
    #include "o8cons.h"

    return;
}