ratfit.c

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

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

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

#define  NDAT  20
#define  ZGR    2
#define  NGR    2
static DOUBLE xdat[NDAT+1],ydat[NDAT+1];
static DOUBLE errbound;
/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    static INTEGER i,j;


    strcpy(name,"ratfitlog");
    n = ZGR+NGR+1;
    nh= 0;
    ng= 2*n+2*NDAT;
    analyt = TRUE ;
    silent = FALSE;
    epsdif = 0.e0;
    del0   = 0.1;
    tau0   = 1.e4;
    errbound=1.0e-1;
    for ( i = 1 ; i<=NDAT ; i++ ) {
      xdat[i] = 0.5+((double) i-1)/((double) 2*NDAT);
      ydat[i] = log(xdat[i]);
    }




    x[1] =  0;
    for ( j=2 ; j <= ZGR+1 ; j++ ) {
      x[j] = -1.0/((double) j-1 );
    }
    for ( j = ZGR+2; j<= n ; j++ ){
      x[j] = 1.0 ;
    }

    for (j = 0 ; j <= 2*NDAT  ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] = 0;
        gunit[3][j] = 0;
    }
    for ( j= 1 ; j<= n; j++ ) {
        gunit[1][j+2*NDAT] = 1 ;
        gunit[2][j+2*NDAT] = j ;
        gunit[3][j+2*NDAT] = 1 ;       /* lower bound */
        gunit[1][j+2*NDAT+n] = 1 ;
        gunit[2][j+2*NDAT+n] = j ;   
        gunit[3][j+2*NDAT+n] = - 1;   /* upper bound */
    }
    return;

}
/***************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    /* change termination criterion */
    te2= TRUE ;  
    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
    icf = icf+1;
    *fx = 0.0 ; 
    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] = 0.0;
    }
    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 j;
    static DOUBLE  numer,denom,argu;

    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;

    if ( i > 2*NDAT ) {
      if ( i > 2*NDAT +n ){
      /*    no upper bound */
        *gxi= 1.0e8-x[i-2*NDAT -n];
        return;
      }
      else {
        j = i-2*NDAT; 
        if ( j > ZGR+1 ) {
          *gxi = x[j] ; 
        }
        else {
          *gxi = x[j]+1.0e8;
        }
        return ;
      }
    }
    else {
    /*   bounds for the approximation error */
      if ( i  <= NDAT ) {
        /*   left bound */
        argu = xdat[i]-1.0;
        denom = x[n]; 
        for ( j= NGR-1; j>=1; j-- ) {
          denom=denom*argu+x[j+ZGR+1] ;
        }
        denom = denom*argu+1.0;
        numer = 0.0 ;
        for ( j= ZGR ; j>=0 ; j-- ) {
          numer = numer*argu+x[j+1];
        }
        *gxi = ydat[i] - numer/denom + errbound;
      }
      else {
        /* right bound */
        denom = x[n];
        argu=xdat[i-NDAT]-1.0;
        for ( j= NGR-1; j>=1; j-- ) { 
          denom=denom*argu+x[j+ZGR+1];
        }
        denom = denom*argu+1.0;
        numer = 0.0 ;
        for ( j= ZGR ; j>=0 ; j-- ) {
          numer = numer*argu+x[j+1];
        }
        *gxi = errbound -ydat[i] + numer/denom ;
      }
      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 gradx[]) {
#define  X extern
#include "o8fuco.h"
#undef   X

    static INTEGER k,j;
    static DOUBLE  numer,denom,argu;
    
    if ( gunit[1][i+nh] != 1 ) cgres[i+nh] = cgres[i+nh]+1;
    if ( i > 2*NDAT ) {
      /* a bound */
      return;
    }
    k =  i ;
    if ( k > NDAT ) {
      k-= NDAT ;
    }
    
    denom = x[n];
    argu = xdat[k]-1.0;
    for ( j= NGR-1; j>=1; j-- ) {   
      denom=denom*argu+x[j+ZGR+1];
    }
    denom = denom*argu+1.0;
    numer = 0.0 ;
    for ( j= ZGR ; j>=0 ; j-- ) {  
      numer = numer*argu+x[j+1];
    }
    if ( i <= NDAT ) { 
      for ( j = 1 ; j<= ZGR+1 ; j++ ) {
         gradx[j] = - pow(argu,j-1)/denom;
      }
      for ( j = 1 ; j<= NGR ; j++ ) {
         gradx[j+ZGR+1] = ((numer/denom)/denom)*pow(argu,j);
      }
    }
    else {
      for ( j = 1 ; j<= ZGR+1 ; j++ ) {
         gradx[j] =  pow(argu,j-1)/denom;
      }
      for ( j = 1 ; j<= NGR ; j++ ) {
         gradx[j+ZGR+1] = -((numer/denom)/denom)*pow(argu,j);
      }
    }  
    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;
}