userfu.c

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

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

DOUBLE A[][10] = 
{
	{0,		0,				0,				0,				0,				0,				0,				0,			0,			  0},
	{0,		0.000374947,	0.000460345,	0.000338638,	-0.000052748,	0.000405596,	0.000437557,	0.000449961, 0.000402754, 0.000460538},
	{0,		0.000000000,	0.000639643,	0.000425192,	0.000078366,	0.000515885,	0.000564882,	0.000580702, 0.000502497, 0.000588405},
	{0,		0.000000000,	0.000000000,	0.000892451,	0.000441222,	0.000391331,	0.000492591,	0.000520105, 0.000421080, 0.000501561},
	{0,		0.000000000,	0.000000000,	0.000000000,	0.003797120,	-0.000099502,	0.000159974,	0.000200895, 0.000081047, 0.000178040},
	{0,		0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000491021,	0.000483901,	0.000493275, 0.000451258, 0.000523641},
	{0,		0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000578161,	0.000583732, 0.000498248, 0.000595467},
	{0,		0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000596655, 0.000507726, 0.000607239},
	{0,		0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000, 0.000454958, 0.000520041},
	{0,		0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000,	0.000000000, 0.000000000, 0.000636609}
};

DOUBLE B3[] = {0,	0.00,	0.00,	0.12,	0.00,	0.00,	0.04,	0.06,	0.00,	0.00};
DOUBLE B5[] = {0,	1.00,	1.00,	0.88,	1.00,	1.00,	0.96,	0.94,	1.00,	1.00};


main() {
#define  X extern
#include "o8comm.h"
#include "o8fint.h"
#undef   X
    void donlp2(void);
    int i;
    DOUBLE sum;
    
    donlp2();
    

    for(i=1; i<=n; ++i)
        printf("%.8f\n", x[i]);

    
	/* B3 should be 0.90006096  */
    for(sum=0, i=1; i<=n; ++i)
        sum += B3[i]*x[i];
    printf("B3 %.8f\n", sum);

	/* B5 should be 0.09993904  */
    for(sum=0, i=1; i<=n; ++i)
        sum += B5[i]*x[i];
    printf("B5 %.8f\n", sum);

    exit(0);
}

#define GENERAL    4	/* this just indicates how many general inequalities
we have. For this example, there's only one general inequality */


void setup0(void) {
#define  X extern
#include "o8comm.h"
#include "o8fint.h"
#undef   X
    
    static INTEGER  i,j;
    
    strcpy(name,"test");


    n  = 9;
    nh = 1;
    ng =   GENERAL    // general inequalities
         + 2*n;        // bounds


    for (i = 1 ; i <= n ; i++) {
        x[i] = 0.e0;
        //x[i] = 1.0/n;
    }


    /*gunit[1][0]=gunit[1][1]=gunit[1][2]=-1;
    gunit[2][0]=gunit[3][0]=gunit[2][1]=gunit[3][1]=gunit[2][2]=gunit[3][2]=0;*/

	for( j=0; j<=(nh+GENERAL); j++)
	{
		gunit[1][j] = -1;
		gunit[2][j] = 0;
		gunit[3][j] = 0;
	}


    for ( j=(nh+GENERAL+1); j<=(nh+ng); j++ ) 
    {
      gunit[1][j]=1; 
      if ( j > (nh+ng-n) )
      {
        gunit[2][j]=j-(nh+ng-n);
        gunit[3][j]=-1;
      }
      else
      {
        gunit[2][j]=j-(nh+GENERAL);
        gunit[3][j]=1;
      }
    }

    tau0=1.e0;		/* how does this (and del0, tau, taubnd, epsdif,
epsfcn) practically affect the solution? */
/* yes, maybe important but not here */

    del0 = 1.e-1;    /* you: 1.0 */
    tau  = 0.1e0;

    taubnd = 1.0;    
    epsdif=1.e-16;
    epsfcn=1.e-16;
    /* modification to use donlp2 high order numerical differentiation : */
    
    analyt   = TRUE;
    difftype = 3;
    
    cold     = TRUE;
    silent   = FALSE;

} 


void setup(void) {
#define  X extern
#include "o8comm.h"
#undef   X
    
    INTEGER i,j;    
    for ( i=1; i<=n; i++)
    { 
        for ( j=i+1; j<=n; j++)
        {
           A[j][i] = A[i][j];
        }
    }
    te2=TRUE;
    te3=TRUE;
    epsx=1.e-7;

    return;
}


/********************************************************/
/*                      objective function              */
/********************************************************/
void ef(DOUBLE x[],DOUBLE *fx) {
#define  X extern
#include "o8fuco.h"
#undef   X

    INTEGER  i, j;
    DOUBLE sum1, sum=0;

    icf = icf+1;

    for(i=1; i<=n; ++i)
    {
        for(j=1, sum1=0; j<=n; ++j)
            sum1 += A[i][j] * x[j];
        sum += x[i] * sum1;
    }
    
    *fx = sum;

    return;
}

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

    INTEGER j, k;
    DOUBLE sum;
    for(j=1; j<=n; ++j)
    {
        for(k=1, sum=0; k<=n; ++k)
            sum += A[j][k] * x[k];
        sum *= 2;

        gradxl[j] = sum;
    }
    
    
    return;
}

DOUBLE sumproduct(DOUBLE x[], DOUBLE B[])
{
#define  X extern
#include "o8fuco.h"
#undef   X
#include "o8cons.h"
    INTEGER j;
    DOUBLE sum;

    for(j=1, sum=0; j<=n; ++j)
        sum += B[j] * x[j];
    return sum;
}

/********************************************************/
/*   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

    int j;
    cres[i] = cres[i]+1;

    /*  sum[ x_i ]  = 1 */
    *hxi = -1;
    for(j=1; j<=n; ++j)
        *hxi += x[j];

    return;
}

/***********************************************************/
/*   compute the gradient of the i-th equality constraint  */
/***********************************************************/
void egradh(INTEGER i,DOUBLE xl[],DOUBLE gradhi[]) {
#define  X extern
#include "o8fuco.h"
#undef   X
#include "o8cons.h"
    INTEGER j;

    for (j = 1 ; j <= n ; j++)
        gradhi[j] = 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
    
    cres[i+nh] = cres[i+nh]+1;

    switch(i) {
    case 1:        *gxi = 1 - sumproduct(x, B5) ;   break;
    case 2:        *gxi = sumproduct(x, B5) ;		break;
    case 3:        *gxi = 1 - sumproduct(x, B3);    break;
    case 4:        *gxi = sumproduct(x, B3) ;		break;

    case 1+GENERAL:    *gxi = x[1] - 0;    break;
    case 2+GENERAL:    *gxi = x[2] - 0;    break;
    case 3+GENERAL:    *gxi = x[3] - 0;    break;
    case 4+GENERAL:    *gxi = x[4] - 0;    break;
    case 5+GENERAL:    *gxi = x[5] - 0;    break;
    case 6+GENERAL:    *gxi = x[6] - 0;    break;
    case 7+GENERAL:    *gxi = x[7] - 0;    break;
    case 8+GENERAL:    *gxi = x[8] - 0;    break;
    case 9+GENERAL:    *gxi = x[9] - 0;    break;

    case 10+GENERAL:    *gxi = 1 - x[1];    break;
    case 11+GENERAL:    *gxi = 1 - x[2];    break;
    case 12+GENERAL:    *gxi = 1 - x[3];    break;
    case 13+GENERAL:    *gxi = 1 - x[4];    break;
    case 14+GENERAL:    *gxi = 1 - x[5];    break;
    case 15+GENERAL:    *gxi = 1 - x[6];    break;
    case 16+GENERAL:    *gxi = 1 - x[7];    break;
    case 17+GENERAL:    *gxi = 1 - x[8];    break;
    case 18+GENERAL:    *gxi = 1 - x[9];    break;
    
    default:    
            break;
    }


    return;
}


void egradg(INTEGER i,DOUBLE x[],DOUBLE gradgi[]) {
#define  X extern
#include "o8fuco.h"
#undef   X
#include "o8cons.h"
    INTEGER j;

    for ( j = 1 ; j <= n ; j++)    
        gradgi[j] = 0;

    switch(i) {
    case 4:    
		for ( j = 1 ; j <= n ; j++)    
			gradgi[j] = B3[j];
		break;
			
    case 2:    
		for ( j = 1 ; j <= n ; j++)    
			gradgi[j] = B5[j];
        break;
    case 3:    
		for ( j = 1 ; j <= n ; j++)    
			gradgi[j] = -B3[j];
        break;

    case 1:    
		for ( j = 1 ; j <= n ; j++)    
			gradgi[j] = -B5[j];
        break;
            
    case 1+GENERAL:        gradgi[1] = 1;        break;
    case 2+GENERAL:        gradgi[2] = 1;        break;
    case 3+GENERAL:        gradgi[3] = 1;        break;
    case 4+GENERAL:        gradgi[4] = 1;        break;
    case 5+GENERAL:        gradgi[5] = 1;        break;
    case 6+GENERAL:        gradgi[6] = 1;        break;
    case 7+GENERAL:        gradgi[7] = 1;        break;
    case 8+GENERAL:        gradgi[8] = 1;        break;
    case 9+GENERAL:        gradgi[9] = 1;        break;

    case 10+GENERAL:    gradgi[1] = -1;        break;
    case 11+GENERAL:    gradgi[2] = -1;        break;
    case 12+GENERAL:    gradgi[3] = -1;        break;
    case 13+GENERAL:    gradgi[4] = -1;        break;
    case 14+GENERAL:    gradgi[5] = -1;        break;
    case 15+GENERAL:    gradgi[6] = -1;        break;
    case 16+GENERAL:    gradgi[7] = -1;        break;
    case 17+GENERAL:    gradgi[8] = -1;        break;
    case 18+GENERAL:    gradgi[9] = -1;        break;
    
    default:    
            break;
    }


    return;

}



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

    return;
}

void solchk(void) {
#define  X extern
#include "o8comm.h"
#undef   X
#include "o8cons.h"

    return;
}