betting.c

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

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

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

/* example Murty 10.15 (Murty K.G.: linear complementarity, linear and nonlinear*/
/* programming, Heldermann, berlin 1988 )                                       */

static DOUBLE bb,qq,pp,pi[9],qi[9];

/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    static INTEGER i,j;
    static DOUBLE  ww,wi[9];

    for (i = 1 ; i <= 8 ; i++) {
        x[i] = i;
    }
    pi[1] =  4.e3;
    pi[2] =  4.e3;
    pi[3] =  4.e3;
    pi[4] =  8.e3;
    pi[5] =  3.e3;
    pi[6] =  8.e3;
    pi[7] = 13.e3;
    pi[8] =  5.e3;
    
    bb   = 5.e2;
    qq   = .83e0;
    
    n    = 8;
    nh   = 0;
    ng   = 9;
    del0 = 1.00e0;
    tau0 = 1.e0;
    tau  = 1e0;
    
    for (j = 0 ; j <= 1 ; j++) {
        gunit[1][j] = -1;
        gunit[2][j] =  0;
        gunit[3][j] =  0;
    }
    for (j = 2 ; j <= 9 ; j++) {
        gunit[1][j] = 1;
        gunit[2][j] = j-1;
        gunit[3][j] = 1;
    }
    wi[1] = 1.e4;
    wi[2] = 1.5e4;
    wi[3] = 5.e3;
    wi[4] = 3.5e4;
    wi[5] = 5.e3;
    wi[6] = 1.e4;
    wi[7] = 1.8e4;
    wi[8] = 1.2e4;

    strcpy(name,"betting8");
     
    silent = FALSE;
    nreset = n;
    ww     = 0.e0;
    pp     = 0.e0;
    analyt = TRUE;
    epsdif = 0.e0;
    for (i = 1 ; i <= n ; i++) {
        ww = ww+wi[i];
        pp = pp+pi[i];
    }
    for (i = 1 ; i <= n ; i++) {
        qi[i] = wi[i]/ww;
    }
    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 INTEGER i,j,l;
    static DOUBLE  sx,xij,fij,term;

    icf  = icf+1;
    term = pp;
    for (i = 1 ; i <= n ; i++) {
        term = term+x[i];
    }
    term = term*qq;
    *fx  = 0.e0;
    for (i = 1 ; i <= n ; i++) {
        for (j = 1 ; j <= n ; j++) {
            if ( j != i ) {
                xij = x[i]/(x[i]+pi[i])+x[j]/(x[j]+pi[j]);
                sx  = 0.e0;
                for (l = 1 ; l <= n ; l++) {
                    if ( l != i && l !=j ) sx = sx+x[l];
                }
                fij = .5e0*(term-x[i]-x[j]-pi[i]-pi[j])*xij+bb-sx;
                *fx = *fx-qi[i]*qi[j]/(1.e0-qi[i])*log(fij);
            }
        }
    }
    return;
}

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

    static INTEGER i,j,k;
    static DOUBLE  sx,fij;
    static DOUBLE  grfij[NX+1],gradt1[NX+1],gradt2[NX+1],grarg[NX+1];
    static DOUBLE  t1,t2,arg;
    
    icgf = icgf+1;
    sx   = 0.e0;
    for (i = 1 ; i <= n ; i++) {
        sx       = sx+x[i];
        gradf[i] = 0.e0;
    }
    for (i = 1 ; i <= n ; i++) {
        for (j = 1 ; j <= n ; j++) {
            if ( i != j ) {
                for (k = 1 ; k <= n ; k++) {
                    grarg[k]  = 0.e0;
                    gradt1[k] = 0.e0;
                    gradt2[k] = 0.e0;
                }
                t2        = x[i]/(x[i]+pi[i])+x[j]/(x[j]+pi[j]);
                gradt2[i] = pi[i]/pow(x[i]+pi[i],2);
                gradt2[j] = pi[j]/pow(x[j]+pi[j],2);
                t1        = (qq*(pp+sx)-x[i]-x[j]-pi[i]-pi[j])*.5e0;
                for (k = 1 ; k <= n ; k++) {
                    gradt1[k] = qq*.5e0;
                }
                gradt1[i] = gradt1[i]-.5e0;
                gradt1[j] = gradt1[j]-.5e0;
                for (k = 1 ; k <= n ; k++) {
                    grfij[k] = t2*gradt1[k]+t1*gradt2[k];
                }
                fij = t1*t2;
                arg = fij+bb-sx+x[i]+x[j];
                for (k = 1 ; k <= n ; k++) {
                    grarg[k] = grfij[k]-1.e0;
                }
                grarg[i] = grarg[i]+1.e0;
                grarg[j] = grarg[j]+1.e0;
                for (k = 1 ; k <= n ; k++) {
                    gradf[k] = gradf[k]-qi[i]*qi[j]/(1.e0-qi[i])/arg*grarg[k];
                }
            }
        }
    }
    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

    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;

    if ( i == 1 ) {
        *gxi = bb;
        for (j = 1 ; j <= n ; j++) {
            *gxi = *gxi-x[j];
        }
    } else {
        *gxi = x[i-1];
    }
    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;

    if ( i == 1 ) {
        for (j = 1 ; j <= n ; j++) {
            gradgi[j] = -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;
}