cuong_nguyen_140503.c

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

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

/* Ignore the array index = 0....only using from index=1 to n to coincide with x indices */
/* ie. ignore A[0][i] and A[i][0]  */

DOUBLE A[4][4] = 
{
    {0,        0,                0,                0,},
    {0,        0.000374947,    0.000460345,    -0.000052748},
    {0,        0.000000000,    0.000639643,    0.000078366},
    {0,        0.000000000,    0.000000000,    0.003797120}
};

/* Again, ignore B1[0], B2[0], etc.  */

            //pad 1st column
DOUBLE B1[] = {0,    0.00,    0.00,    0.00};
DOUBLE B2[] = {0,    0.00,    0.00,    0.00};
DOUBLE B3[] = {0,    0.35,    0.00,    0.00};
DOUBLE B4[] = {0,    0.00,    0.00,    0.00};
DOUBLE B5[] = {0,    0.65,    1.00,    1.00};
DOUBLE B6[] = {0,    0.00,    0.00,    0.00};
DOUBLE B7[] = {0,    0.00,    0.00,    0.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]);

    
    for(sum=0, i=1; i<=n; ++i)
        sum += B2[i]*x[i];
    printf("B2 %.8f\n", sum);
/* CHANGED */
/*  since general =1 it makes no sense to test two inequalities here
    for(sum=0, i=1; i<=n; ++i)
        sum += B5[i]*x[i];
    printf("B5 %.8f\n", sum);
*/
    exit(0);
}

/*
****************************************************************************
*/
/*                              donlp2 standard setup
*/
/*
****************************************************************************
*/

#define GENERAL    1	/* 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  = 3;
    nh = 1;
    ng =   GENERAL    // general inequalities
         + 2*n;        // bounds


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

/*
#define G(i, x, y, z)    {gunit[1][i]=x; gunit[2][i]=y; gunit[3][i]=z;}
    
    G(0, -1, 0, 0);            // fx
    G(1, -1, 0, 0);            // first equality (sum(x) = 1)

    G(2, -1, 0, 0);            // general inequality


    G(2+GENERAL, 1, 1, 1);        // lower bounds
    G(3+GENERAL, 1, 2, 1);        // lower bounds
    G(4+GENERAL, 1, 3, 1);        // lower bounds

    G(5+GENERAL, 1, 1, -1);    // upper bounds
    G(6+GENERAL, 1, 2, -1);
    G(7+GENERAL, 1, 3, -1);
*/
    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=3; j<=8; j++ ) 
    {
      gunit[1][j]=1; 
      if ( j > 5 )
      {
        gunit[2][j]=j-5;
        gunit[3][j]=-1;
      }
      else
      {
        gunit[2][j]=j-2;
        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.e0;    /* you: 0.1 */
    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;

} 


/*
****************************************************************************
*/
/*                                 special setup
*/
/*
****************************************************************************
*/
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;
    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

    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;
}

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

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

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

void gsumit(DOUBLE grad[], DOUBLE A[], DOUBLE factor)
{
#define  X extern
#include "o8fuco.h"
#undef   X
#include "o8cons.h"
    INTEGER j;

    for(j=1; j<=n; ++j)
        grad[j] = factor * A[j-1];
}

/*
****************************************************************************
*/
/*              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 =sumproduct(x, B2);    break;
/*    B2 was zero vector and you wanted B3 (and B5)*/
/*    CHANGED FROM B2 TO B3 */    
    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 = 1 - x[1];    break;
    case 5+GENERAL:    *gxi = 1 - x[2];    break;
    case 6+GENERAL:    *gxi = 1 - x[3];    break;
    
    default:    
            break;
    }


    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;

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

    switch(i) {
/* CHANGED FROM B2 TO B3 */
    case 1:    
            gradgi[1] = B2[1];
            gradgi[2] = B2[2];
            gradgi[3] = B2[3];
            break;
            
    case 1+GENERAL:        gradgi[1] = 1;        break;
    case 2+GENERAL:        gradgi[2] = 1;        break;
    case 3+GENERAL:        gradgi[3] = 1;        break;

    case 5+GENERAL:    gradgi[1] = -1;        break;
    case 6+GENERAL:    gradgi[2] = -1;        break;
    case 7+GENERAL:    gradgi[3] = -1;        break;
    
    default:    
            break;
    }


    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;
}