parametric.c

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

/* **************************************************************************** */
/*                                 user functions                               */
/* **************************************************************************** */
#include "o8para.h"
/* **************************************************************************** */
/* this is a modified version of a problem                                      */
/* supplied by Feng Ling                                                        */
/* we show here how to use donlp2 in a parametric application:                  */
/* the bounds of constraining box are thightened until infeasibility            */
/* of the problem is obtained                                                   */
/* in main a factor bdfac is defined which enters into the definition           */
/* of the bounds                                                                */
/* **************************************************************************** */
static DOUBLE   alamda[4],f0[4],c[4],
                fmag,fdelta[4],fdeltr[1],
                fdeltb[4],dm0,wrench,xt[4],yt[4],
                xbal,ybal,bdfac;
static INTEGER  case_;

main() {
    #define  X extern
    #include "o8comm.h"
    #undef   X

    void donlp2(void);

    case_ = 0;
    bdfac = 50.e0;

    while ( bdfac > 1.0e0 ) {
        case_ = case_+1;
        bdfac = max(1.e0, bdfac-1.e0);

        donlp2();
    }
    exit(0);
}

/* **************************************************************************** */
/*                              donlp2 standard setup                           */
/* **************************************************************************** */
void setup0(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    
    void setpar(void);
    
    static INTEGER  j,nuser = 6;
    static DOUBLE   xst0[7];
    static LOGICAL  first = TRUE;

    /* case_ is increased by main. Any case_ has its own output-file */

    sprintf(name,"forces%2.2i",case_);
    
    /* x is initial guess and also holds the current solution */
    /* problem dimension n = dim(x), nh = dim(h), ng = dim(g) */
    
    n      = nuser;
    nh     = 3;
    ng     = 9;
    analyt = TRUE;
    epsdif = 0.e0;

    /* del0 and tau0: see below */
    
    del0 = 1.e0;
    tau0 = 1.e-1;
    if ( first ) {
        cold = TRUE;
        
        setpar();
        
        for (j = 1 ; j <= n ; j++) {
            if (j >3) xst0[j] = alamda[j-3];
            if (j<=3) xst0[j] = 500.e0;
            x[j] = xst0[j];
        }
        for (j = 0 ; j <= 3 ; j++) {
            gunit[1][j] = -1;
            gunit[2][j] = 0;
            gunit[3][j] = 0;
        }
        for (j = 4 ; j <= 9 ; j++) {
            gunit[1][j] = 1;
            gunit[2][j] = j-3;
            gunit[3][j] = 1;
        }
        for (j = 10 ; j <= 12 ; j++) {
            gunit[1][j] = 1;
            gunit[2][j] = j-6;
            gunit[3][j] = -1;
        }
        first = FALSE;
    } else {
        itstep = 0;
        cold   = FALSE;
    }
    return;
}

/* **************************************************************************** */
/*                                 special setup                                */
/* **************************************************************************** */
void setup(void) {
    #define  X extern
    #include "o8comm.h"
    #undef   X
    /* te0=TRUE;
    /* giving short protocol on stdout                                  */
    /* any change of donlp2's standard parameter initialization must be */
    /* done here                                                        */

    return;
}

void setpar(void) {

    /* special user data entering the function definiton */
    
    #define  X extern
    #include "o8comm.h"
    #undef   X

    static INTEGER j;

    for (j = 1 ; j <= 3 ; j++) {
        alamda[j] = 0.01e0;
        c[j]      = -1.e1*j;
        fdelta[j] = 0.1e0;
        fdeltb[j] = 0.05e0;
    }
    xbal      = 0.e0;
    ybal      = 0.e0;
    xt[1]     = -1.e0;
    yt[1]     = 1.e0;
    xt[2]     = 0.e0;
    yt[2]     = -1.e0;
    xt[3]     = 1.e0;
    yt[3]     = 1.e0;
    f0[1]     = 10.e0;
    f0[2]     = 20.e0;
    f0[3]     = 5.e0;
    fmag      = 1.e-1;
    dm0       = 2.e0;
    wrench    = 1.e0;
    fdeltr[0] = 0.5e0;
    
    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 j;
    
    icf = icf+1;
    *fx = 0.e0;
    for (j = 1 ; j <= 3 ; j++) {
        *fx = *fx+pow(f0[j]+x[j],2)*(1+pow(x[j+3],2)*pow(c[j],2));
    }
    return;
}

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

    static INTEGER j;
    
    icgf = icgf+1;
    for (j = 1 ; j <= 3 ; j++) {
        gradf[j]   = 2.e0*(1.e0+pow(x[j+3],2)*pow(c[j],2))*(f0[j]+x[j]);
        gradf[j+3] = 2.e0*pow(c[j],2)*pow(f0[j]+x[j],2)*x[j+3];
    }
    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

    static DOUBLE   ati[4],pc,ps,pss1,pss2,pss3,pss4;
    static INTEGER  j;
    
    cres[i] = cres[i]+1;
    *hxi    = 0.e0;
    switch (i) {
    case 1:
        pc = fmag*cos(fdeltr[0]);
        for (j = 1 ; j <= 3 ; j++) {
            ati[j] = cos(fdelta[j])-x[j+3]*c[j]*sin(fdelta[j]);
        }
        for (j = 1 ; j <= 3 ; j++) {
            *hxi = *hxi+ati[j]*(x[j]+f0[j]);
        }
        *hxi = *hxi+pc;
        
        break;
    case 2:
        ps = fmag*sin(fdeltr[0]);
        for (j = 1 ; j <= 3 ; j++) {
            ati[j] = sin(fdelta[j])+x[j+3]*c[j]*cos(fdelta[j]);
        }
        for (j = 1 ; j <= 3 ; j++) {
            *hxi = *hxi+ati[j]*(x[j]+f0[j]);
        }
        *hxi = *hxi+ps;
        
        break;
    case 3:
        for (j = 1 ; j <= 3 ; j++) {
            pss1 = sqrt(1+pow(c[j],2)*pow(x[j+3],2));
            pss2 = f0[j]+x[j];
            pss3 = sqrt(pow(ybal-yt[j],2)+pow(xbal-xt[j],2));
            pss4 = sin(fdelta[j]-fdeltb[j]+atan(x[j+3]*c[j]));
            *hxi = *hxi+pss1*pss2*pss3*pss4;
        }
        *hxi = *hxi-dm0*wrench;
        
        break;
    }
    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

    static DOUBLE   pss1,pss2,pss3,pss4;
    static DOUBLE   gpss1[7],gpss2[7],gpss4[7],arg;
    static INTEGER  j,k;
    
    if ( gunit[1][i] != 1 ) cgres[i] = cgres[i]+1;

    switch (i) {
    case 1:
        for (j = 1 ; j <= 3 ; j++) {
            gradhi[j]   = cos(fdelta[j])-x[j+3]*c[j]*sin(fdelta[j]);
            gradhi[j+3] = -c[j]*sin(fdelta[j])*(x[j]+f0[j]);
        }
        break;
    case 2:
        for (j = 1 ; j <= 3 ; j++) {
            gradhi[j]   = sin(fdelta[j])+x[j+3]*c[j]*cos(fdelta[j]);
            gradhi[j+3] = c[j]*cos(fdelta[j])*(x[j]+f0[j]);
        }
        break;
    case 3:
        for (k = 1 ; k <= 6 ; k++) {
            gradhi[k] = 0.e0;
        }
        for (j = 1 ; j <= 3 ; j++) {
            pss1 = sqrt(1+pow(c[j],2)*pow(x[j+3],2));
            pss2 = f0[j]+x[j];
            pss3 = sqrt(pow(ybal-yt[j],2)+pow(xbal-xt[j],2));
            arg  = fdelta[j]-fdeltb[j]+atan(x[j+3]*c[j]);
            pss4 = sin(arg);
            for (k = 1 ; k <= 6 ; k++) {
                gpss1[k] = 0.e0;
                gpss2[k] = 0.e0;
                gpss4[k] = 0.e0;
            }
            gpss1[j+3] = x[j+3]*pow(c[j],2)/pss1;
            gpss2[j]   = 1.e0;
            gpss4[j+3] = cos(arg)*c[j]/(1.e0+pow(x[j+3]*c[j],2));
            for (k = 1 ; k <= 6 ; k++) {
                gradhi[k] = gradhi[k]+pss3*(gpss1[k]*pss2*pss4+
                            pss1*gpss2[k]*pss4+
                            pss1*pss2*gpss4[k]);
            }
        }
        break;
    }
    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 k;
        
    if ( gunit[1][i+nh] == -1 ) cres[i+nh] = cres[i+nh]+1;
    k = (i-1)/3;
    k = k+1;
    switch (k) {
    case 1:
        *gxi = x[i]+f0[i]*bdfac;
        break;
    case 2:
        *gxi = x[i]+bdfac;
        break;
    case 3:
        *gxi = bdfac-x[i-3];
        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

    static INTEGER k,j;
    
    if ( gunit[1][i+nh] != 1 ) cgres[i+nh] = cgres[i+nh]+1;
    for (j = 1 ; j <= 6 ; j++) {
        gradgi[j] = 0.e0;
    }
    k = (i-1)/3;
    k = k+1;
    switch (k) {
    case 1:
    case 2:
        gradgi[i] = 1.e0;
        break;
    case 3:
        gradgi[i-3] = -1.e0;
        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;
}