iwlaustral_exa.c

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

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

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

/*
****************************************************************************
*/
/*                              donlp2 standard setup
*/
/*
****************************************************************************
*/
void setup0(void) {
#define  X extern
#include "o8comm.h"
#include "o8fint.h"
#undef   X
    
    static INTEGER  i;
    
    /* name is ident of the example/user and can be set at users will
*/
    /* the first static character must be alphabetic. 40 characters maximum
*/

    strcpy(name,"test");

    n  = 7;
    nh = 2;
    ng = 15;
    for (i = 1 ; i <= n ; i++) {
        //x[i] = 0.e0;
		x[i] = 1.0/n;
    }
    for ( i=0; i<=3; i++ )
    { /* two equality constraints and one general inequality*/
        gunit[1][i] = -1;
        gunit[2][i] = 0; 
        gunit[3][i] = 0; 
    }

    for ( i=4; i<=17; i++ )
    { /* the fourteen bounds in [0,1]*/
        gunit[1][i] = 1 ;
        if ( i <= 10 )
        {
          gunit[2][i] = i-3;
          gunit[3][i] = 1;
        }
        else
        {
          gunit[2][i] = i-10;
          gunit[3][i] = -1;
        }
     }

    tau0=1.e0;

    del0 = 1.e-1;
    tau  = 0.1e0;

    taubnd = 1.0;    
    epsdif=1.e-16;
    epsfcn=1.e-16;
    
    analyt   = TRUE;
    difftype = 3;
    
    cold     = TRUE;
    silent   = FALSE;

    return;
}

/***************************************************************************/
/*                                 special setup                           */
/***************************************************************************/



/* NEW */
DOUBLE amat[8][8], emat[8][8];

DOUBLE fOPRCoVar[][8] = 
{
	{0.00000000,	0.000000000,	0.000000000,	0.000000000,
0.000000000,	0.000000000,	0.000000000,	0.000000000},
	{0.00000000,	0.002962261,	0.002648278,	0.002275971,
0.001138859,	0.001764342,	0.000345131,	-0.000413952},
	{0.00000000,	0.00000000,		0.007321681,	0.002586208,
0.000258487,	0.001732468,	0.001955888,	0.000478739},
	{0.00000000,	0.00000000,		0.000000000,	0.003586436,
0.000929674,	0.002487289,	0.001062284,	0.000251558},
	{0.00000000,	0.00000000,		0.000000000,	0.000000000,
0.004545468,	0.001801269,	0.000584452,	0.000972678},
	{0.00000000,	0.00000000,		0.000000000,	0.000000000,
0.000000000,	0.004790799,	0.002206663,	0.001600869},
	{0.00000000,	0.00000000,		0.000000000,	0.000000000,
0.000000000,	0.000000000,	0.003337734,	0.001533548}
};
DOUBLE fOPR[] = {	0, 0.054427,	0.085567,	0.059887,
0.067420,	0.069216,	0.057773,	0.052981};

DOUBLE fTECoVar[][8] = 
{
	{0.00000000,	0.00000000,	0.000000000,	0.000000000,
0.000000000,	0.000000000,	0.000000000,	0.000000000},
	{0.00000000,	12.5791242,	18.89880695,	14.09595257,
3.219307111,	10.63632112,	2.393812436,	-1.523227897},
	{0.00000000,	0.00000000,	87.80622602,	26.91750756,
1.227931082,	17.55163566,	22.79785214,	2.960452545},
	{0.00000000,	0.00000000,	0.000000000,	32.39600395,
3.832866176,	21.86933125,	10.74601498,	1.350064545},
	{0.00000000,	0.00000000,	0.000000000,	0.000000000,
8.553351753,	7.228573138,	2.698484422,	2.382592159},
	{0.00000000,	0.00000000,	0.000000000,	0.000000000,
0.000000000,	41.00135368,	21.72819457,	8.362817803},
	{0.00000000,	0.00000000,	0.000000000,	0.000000000,
0.000000000,	0.000000000,	37.81263163,	9.217045062}
};
DOUBLE fTE[] =  {	0, 3.546706,	9.370498,	5.691749,
2.924611,	6.403230,	6.149198,	2.991706};





void setup(void) {
#define  X extern
#include "o8comm.h"
#undef   X
INTEGER i,j;    
    /*epsx = 1.e-8;*/
    for ( i=1; i<=n; i++)
      { amat[i][i] = pow(fOPR[i],2);
        emat[i][i] = pow(fTE[i],2);
        for ( j=i+1; j<=n; j++)
        {
           amat[i][j] = amat[j][i] = fOPRCoVar[i][j];
           emat[i][j] = emat[j][i] = fTECoVar[i][j]; 
        }
      }

    return;
}

/*************************************************************************/
/*  the user may add additional computations using the computed solution */
/*************************************************************************/
void solchk(void) {
#define  X extern
#include "o8comm.h"
#undef   X
#include "o8cons.h"

    return;
}

/**************************************************************************/
/*                               objective function                       */
/**************************************************************************/
DOUBLE fError[] = {	0, 0.177088,	0.056855,	0.181746,
0.112348,	0.127406,	0.118753,	0.162796};

/**************************************************************************/
/*                      objective function                                */
/***************************************************************************/ 


void ef(DOUBLE x[],DOUBLE *fx) {
#define  X extern
#include "o8fuco.h"
#undef   X
    DOUBLE sum,sum1;
    INTEGER i,j;
    icf += 1;

    /* maximize quadratic function x'*amat*x  */
    sum=0.0;
    for  (i=1; i<=n; i++)
    {
      sum1=0.0;
      for (j=1; j<=n; j++)
      {
         sum1 += amat[i][j]*x[j];
      }
      sum += x[i]*sum1;
    }
    *fx = -sum;
    return;
}

/***************************************************************************/
/*                          gradient of objective function                 */
/***************************************************************************/
void egradf(DOUBLE x[], DOUBLE gradf[]) {
#define  X extern
#include "o8fuco.h"
#undef   X
#include "o8cons.h"
        DOUBLE sum;
	INTEGER  i,j;
	for ( i=1 ; i<=n ; i++ )
        {
		sum=0.0;
                for ( j=1; j<=n; j++)
                {
                   sum -= amat[i][j]*x[j];
                }
                gradf[i] = 2.0*sum;
        }
                

    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

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

	switch(i) {
	case 1:
		/*  sum[ x_i ]  = 1 */
		*hxi = -1.0;
		for(j=1; j<=n; ++j)
			*hxi += x[j];
		break;
	
	case 2:
		/*  sum[ x_i * fError_i * 100 ]  = 10  */
		*hxi = -10;
		for(j=1; j<=n; ++j)
			*hxi += x[j] * fError[j] * 100;
		break;

	default:
		fprintf(stderr, "undefined case in eh\n");
                exit(1);
	}
    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;

    switch(i) {
    case 1:
		for (j = 1 ; j <= n ; j++)
			gradhi[j] = 1;
        break;
    case 2:
		for (j = 1 ; j <= n ; j++)
			gradhi[j] = fError[j]*100;
        break;
	default:
		fprintf(stderr, "undefined case in  egradh\n");
                exit(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
        INTEGER j,k;
        DOUBLE sum1,sum2,sum3;

/*  	THE NEW VERSION */
/*      i=1 represents the general constraint, i=2 to 2*n+1 the 2n lower
        and upper bounds */
        if ( i == 1 ) 
        {
          cres[i+nh] += 1;
          sum1=11.45985986;
/*  BETTER USE A NAMED PARAMETER HERE, WHICH MAY BE SET OR
    LATER ON READ IN FROM A FILE IN SETUP0 
*/
          
          sum2=0.0;
          for ( j=1; j<=n; j++)
          {
            sum1 -= 100.0*fError[j]*x[j];
            sum3=0.0;
            for (k=1; k<=n; k++)
            {  
              sum3 += emat[j][k]*x[k];
            }
            sum2 += x[j]*sum3;
/* sum2 is the quadratic form built from emat */
          }
          *gxi = sum2+sum1;
   
          return;
        }
        else
        {
           j=i-1;
           if ( j <= n )
           {
              *gxi = x[j];
              return;
           }
           else
           {
              *gxi=1.0-x[j-n];
              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
#include "o8cons.h"
	INTEGER j,k;
        DOUBLE sum;

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

     if ( i==1 ) 
     {
       cgres[i+nh] += 1;
       for ( j=1; j<=n; j++ )
       { 
         sum=0.0;
         for ( k=1; k<=n; k++)
         {
           sum += emat[j][k]*x[k];
         }
         gradgi[j] = 2.0*sum-100.0*fError[j];
       }
       return;
     }
     else
     {
        j=i-1;
        if ( j<=n )
        {
          gradgi[j] = 1.0;
          return;
        }
        else
        {
          gradgi[j] = -1.0;
          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;
}