sqpsolve.src

没有说明

/*
** sqpSolve.src
**
**
** (C) Copyright 1997  Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
******************************************************************************
**
**> sqpSolve
**
**  Purpose:  solve the nonlinear programming problem
**
**  Format:   { x,f,lagr,retcode } = sqpSolve(&fct,start)
**
**
**  Input:    &fct      pointer to a procedure that computes the function to
**                      be minimized.  This procedure must have one input
**                      argument, a vector of parameter values, and one
**                      output argument, the value of the function evaluated
**                      at the input vector of parameter values.
**
**           start      Kx1 vector of start values
**
**
**  Output:      x      Kx1 vector of parameters at minimum
**
**               f      scalar, function evaluated at x
**
**            lagr      vector, created using VPUT.  Contains the Lagrangean
**                      for the constraints.  The may be extracted with the
**                      VREAD command using the following strings:
**
**                        "lineq"    - Lagrangeans of linear equality
                                          constraints,
**                        "nlineq"   - Lagrangeans of nonlinear equality
                                          constraints
**                        "linineq"  - Lagrangeans of linear inequality
                                          constraints
**                        "nlinineq" - Lagrangeans of nonlinear inequality
**                                     constraints
**                        "bounds"   - Lagrangeans of bounds
**
**                      Whenever a constraint is active, its associated
**                      Lagrangean will be nonzero.
**
**
**         retcode      return code:
**
**                         0   normal convergence
**                         1   forced exit
**                         2   maximum number of iterations exceeded
**                         3   function calculation failed
**                         4   gradient calculation failed
**                         5   Hessian calculation failed
**                         6   line search failed
**                         7   error with constraints
**
**  Input Globals:
**
**   _sqp_A            MxK matrix, linear equality constraint coefficients
**   _sqp_B            Mx1 vector, linear equality constraint constants
**
**                      These globals are used to specify linear equality
**                      constraints of the following type:
**
**                         _sqp_A * X = _sqp_B
**
**                      where X is the Kx1 unknown parameter vector.
**
**   _sqp_EqProc       scalar, pointer to a procedure that computes
**                      the nonlinear equality constraints.  For example,
**                      the statement:
**
**                          _sqp_EqProc = &eqproc;
**
**                      tells CO that nonlinear equality constraints
**                      are to be placed on the parameters and where the
**                      procedure computing them is to be found.
**                         The procedure must have one input argument, the
**                      Kx1 vector of parameters, and one output argument,
**                      the Rx1 vector of computed constraints that are
**                      to be equal to zero.  For example, suppose that
**                      you wish to place the following constraint:
**
**                         P[1] * P[2] = P[3]
**
**                      The proc for this is:
**
**                         proc eqproc(p);
**                            retp(p[1]*[2]-p[3]);
**                         endp;
**
**   _sqp_C             MxK matrix, linear inequality constraint coefficients
**   _sqp_D             Mx1 vector, linear inequality constraint constants
**
**                       These globals are used to specify linear inequality
**                       constraints of the following type:
**
**                         _sqp_C * X >= _sqp_D
**
**                       where X is the Kx1 unknown parameter vector.
**
**   _sqp_IneqProc     scalar, pointer to a procedure that computes
**                      the nonlinear inequality constraints.  For example
**                      the statement:
**
**                           _sqp_EqProc = &ineqproc;
**
**                      tells CO that nonlinear equality constraints
**                      are to be placed on the parameters and where the
**                      procedure computing them is to be found.
**                      The procedure must have one input argument, the
**                      Kx1 vector of parameters, and one output argument,
**                      the Rx1 vector of computed constraints that are
**                      to be equal to zero.  For example, suppose that
**                      you wish to place the following constraint:
**
**                         P[1] * P[2] >= P[3]
**
**                      The proc for this is:
**
**                         proc ineqproc(p);
**                            retp(p[1]*[2]-p[3]);
**                         endp;
**
**   _sqp_Bounds       Kx2 matrix, bounds on parameters.  The first column
**                      contains the lower bounds, and the second column the
**                      upper bounds.  If the bounds for all the coefficients
**                      are the same, a 1x2 matrix may be used.
**                      Default = { -1e256 1e256 }
**
**   _sqp_GradProc     scalar, pointer to a procedure that computes the
**                      gradient of the function with respect to the
**                      parameters.  For example, the statement:
**
**                            _sqp_GradProc=&gradproc;
**
**                      tells CO that a gradient procedure exists as well
**                      where to find it.  The user-provided procedure has
**                      two input arguments, a Kx1 vector of parameter values
**                      and an NxP matrix of data.  The procedure returns a
**                      single output argument, an NxK matrix of gradients
**                      of the log-likelihood function with respect to the
**                      parameters evaluated at the vector of parameter values.
**
**                      Default = 0, i.e., no gradient procedure has been
**                      provided.
**
**   _sqp_HessProc     scalar, pointer to a procedure that computes the
**                      hessian, i.e., the matrix of second order partial
**                      derivatives of the function with respect to the
**                      parameters.  For example, the instruction:
**
**                           _sqp_HessProc=&hessproc;
**
**                      will tell OPTMUM that a procedure has been provided
**                      for the computation of the hessian and where to find
**                      it.  The procedure that is provided by the user must
**                      have two input arguments, a Px1 vector of parameter
**                      values and an NxK data matrix.  The procedure returns
**                      a single output argument, the PxP symmetric matrix of
**                      second order derivatives of the function evaluated at
**                      the parameter values.
**
**
**
**   _sqp_MaxIters     scalar, maximum number of iterations. Default = 1e+5.
**                     Termination can be forced by pressing C on the keyboard
**
**   _sqp_DirTol       scalar, convergence tolerance for gradient of estimated
**                      coefficients.  Default = 1e-5.  When this criterion has
**                      been satisifed NLPSolve will exit the iterations.
**
**   _sqp_ParNames     Kx1 character vector,  parameter names
**
**   _sqp_PrintIters   scalar, if nonzero, prints iteration information.
**                     Default = 0.  Can be toggled during iterations by
**                     pressing P on the keyboard.
**
**   _sqp_FeasibleTest scalar, if nonzero, parameters are tested for feasibility
**                      before computing function in line search.  If function
**                      is defined outside inequality boundaries, then this test
**                      can be turned off.
**
**   _sqp_RandRadius   scalar, If zero, no random search is attempted.  If
**                     nonzero, it is the radius of random search which is
**                     invoked whenever the usual line search fails.
**                     Default = .01.
**
**    __output          scalar, if nonzero, results are printed.  Default = 0.
**
**
**  Remarks:  Pressing C on the keyboard will terminate iterations, and
**            pressing P will toggle iteration output.
**
**            NLPSolve is recursive, that is, it can call a version of itself
**            with another function and set of global variables,
*/

#include gauss.ext
#include sqpsolve.ext



proc (4) = sqpSolve(fnct,start);
    local x,f,retcode,lagrange,iters,t0,t1,lbl,fmt;

    if __output;
        t0 = hsec;
    endif;
    { x,f,lagrange,retcode,iters } = _sqpsolve(fnct, start,
        _sqp_GradProc, _sqp_HessProc, _sqp_DirTol,
        _sqp_MaxIters, _sqp_PrintIters, _sqp_FeasibleTest,
        _sqp_A, _sqp_B, _sqp_C, _sqp_D, _sqp_Bounds, _sqp_EqProc,
        _sqp_IneqProc, _sqp_RandRadius );


    if __output;

        t1 = hsec;
        print;
        call header("NLPSolve","",_rtl_ver);
        print;
        print "return code = " ftos(retcode,"%*.*lf",4,0);
        if retcode == 0;
           print "normal convergence";
        elseif retcode == 1;
           print "forced termination";
        elseif retcode == 2;
           print "maximum number of iterations exceeded";
        elseif retcode == 3;
           print "function calculation failed";
        elseif retcode == 4;
           print "gradient calculation failed";
        elseif retcode == 5;
           print "Hessian calculation failed";
        elseif retcode == 6;
           print "line search failed";
        elseif retcode == 7;
           print "error with constraints";
        endif;
        print;
        print "Value of objective function " ftos(f,"%*.*lf",15,6);
        print;
        print "Parameters    Estimates";
        print "-----------------------------------------";

        if _qn_ParNames $== 0;
            lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
        else;
            lbl = _qn_ParNames;
        endif;
        let fmt[2,3] = "-*.*s" 9 8 "*.*lf" 14 4;
        call printfm(lbl~x,0~1,fmt);

        print;
        print;
        print "Linear Equality Lagrangean Coefficients";
        print vread(lagrange,"lineq");
        print;
        print "Linear Inequality Lagrangean Coefficients";
        print vread(lagrange,"linineq");
        print;
        print "Nonlinear Equality Lagrangean Coefficients";
        print vread(lagrange,"nlineq");
        print;
        print "Nonlinear Inequality Lagrangean Coefficients";
        print vread(lagrange,"nlinineq");
        print;
        print "Bounds Lagrangean Coefficients";
        print vread(lagrange,"bounds");

        print;
        print "Number of iterations    " ftos(iters,"%*.*lf",5,0);
        print "Minutes to convergence  " ftos((t1-t0)/6000,"%*.*lf",10,5);
    endif;

    retp(x,f,lagrange,retcode);
endp;



proc (5) = _sqpsolve(&fnct,start,LNSgdprc,LNShsprc,LNSDtol,LNSmiter,
  LNSPiters,LNSFtest,LNS_A,LNS_B,LNS_C,LNS_D,LNS_Bnds,
  LNSeproc,LNSiproc,LNSrteps);

    /* ------- LOCALS ----------- */
    local g,s,h,ky,old,f0,x0,iter,
        np,gproc,hsproc,LNSlagr,k1,numeq,qp_m,
        qp_t,qp_xl,qp_xu,qp_ret,qp_maxit,lagr1,lagr2,qp_a,qp_b,qp_d,qp_lql,
        numNlinEqC,numNlinIneqC,eqproc,ineqproc,fnct:proc,rteps,delta,
        ub,lb,cdelta,dg,tt,f1,f2,w1,sprev,s2prev,rprev,r2prev,j,
        sprev2,s2prev2,sp2,dsprev,vv,zz,ab,a,b,qv,MaxTry;


    clear numEq,lagr1,lagr2,qp_m,qp_ret,LNSlagr;
    clear numNlinEqC,numNlinIneqC;

    MaxTry = 100;
    vv = ones(2,2);
    vv[1,2] = -1;
    zz = zeros(2,1);
    sp2 = zeros(1,2);
    start = vec(start);
    old = ndpcntrl(0,0);
    call ndpcntrl(1,1);

    qp_maxit = 1000;
    qp_d = .01*start;
    qp_t = 1e+256*ones(rows(start),1);
    qp_lql = 1;

    LNSlagr = vput(LNSlagr,error(0),"lineq");
    LNSlagr = vput(LNSlagr,error(0),"linineq");
    LNSlagr = vput(LNSlagr,error(0),"nlineq");
    LNSlagr = vput(LNSlagr,error(0),"nlinineq");
    LNSlagr = vput(LNSlagr,error(0),"bounds");


    if LNSgdprc /= 0;
        gproc = LNSgdprc;
        local gproc:proc;
    endif;

    if LNShsprc /= 0;
        hsproc = LNShsprc;
        local hsproc:proc;
    endif;

    if not scalmiss(LNS_A);
        qp_m = qp_m + rows(LNS_A);
        numEq = numEq + rows(LNS_A);
    endif;

    if not scalmiss(LNS_C);
        qp_m = qp_m + rows(LNS_C);
    endif;

    if not scalmiss(LNSeproc);
        eqproc = LNSeproc;
        local eqproc:proc;
        numNlinEqC = rows(EqProc(start));
        qp_m = qp_m + numNlinEqC;
        numEq = numEq + numNlinEqC;
    endif;

    if not scalmiss(LNSiproc);
        ineqproc = LNSiproc;
        local ineqproc:proc;
        numNlinIneqC = rows(IneqProc(start));
        qp_m = qp_m + numNlinIneqC;
    endif;

    if not scalmiss(LNS_Bnds);
        if cols(LNS_Bnds) /= 2 or (rows(LNS_Bnds) /= rows(start) and
            rows(LNS_Bnds) /= 1);
            if not trapchk(4);
                errorlog "ERROR:  LNS_Bnds is not correctly defined";
            endif;
            call ndpcntrl(old,0xffff);
            ndpclex;
            retp(start,error(0),LNSlagr,7,0);
        endif;
    endif;

    if qp_m == 0;           /* if no constraints other than bounds  */
        qp_m = 1;           /* set to one dummy constraint */
    endif;

    x0 = start;
    qp_d = .1 * x0;
    np = rows(x0);          /* Number of parameters to estimate */


#ifUNIX
    for iter(1,LNSmiter,1);
#else
    iter = 1; do until iter > LNSmiter;
#endif

    f0 = fnct(x0);
    if scalmiss(f0) or (f0 $== __INFp) or (f0 $== __INFn) or (f0 $==
         __INDEFp) or (f0 $== __INDEFn);
            retp(x0,f0,LNSlagr,3,iter);
    endif;


    if LNSgdprc /= 0;
        g = gproc(x0);
    else;
        g = gradp(&fnct,x0)';
    endif;