BCNLS

Solves a nonlinear least-squares problem subject to bounds on the variables and general linear constraints.

Required Arguments

FCN — User-supplied subroutine to evaluate the function to be minimized. The usage is CALL FCN (M, N, X, F), where
M - Number of functions.   (Input)
N - Number of variables.   (Input)
X - Array of length N containing the point at which the function will be evaluated.   (Input)
F - Array of length M containing the computed function at the point X.   (Output)
The routine FCN must be declared EXTERNAL in the calling program.

M — Number of functions. (Input)

C — MCON ´ N matrix containing the coefficients of the MCON general linear constraints. (Input)

BL — Vector of length MCON containing the lower limit of the general constraints. (Input).

BU — Vector of length MCON containing the upper limit of the general constraints. (Input).

IRTYPE — Vector of length MCON indicating the types of general constraints in the matrix C. (Input)
Let R(I) = C(I, 1)*X(1) + … + C(I, N)*X(N). Then the value of IRTYPE(I) signifies the following:

     IRTYPE(I)                         I-th CONSTRAINT
   0                  BL(I).EQ.R(I).EQ.BU(I)
   1                      R(I).LE.BU(I)
   2                               R(I).GE.BL(I)
   3                               BL(I).LE.R(I).LE.BU(I)

XLB — Vector of length N containing the lower bounds on variables; if there is no lower bound on a variable, then 1.0E30 should be set as the lower bound. (Input)

XUB — Vector of length N containing the upper bounds on variables; if there is no upper bound on a variable, then -1.0E30 should be set as the upper bound. (Input)

X — Vector of length N containing the approximate solution. (Output)

Optional Arguments

N — Number of variables. (Input)
Default: N = size (C,2).

MCON — The number of general linear constraints for the system, not including simple bounds. (Input)
Default: MCON = size (C,1).

LDC — Leading dimension of C exactly as specified in the dimension statement of the calling program. (Input)
LDC must be at least MCON.
Default: LDC = size (C,1).

XGUESS — Vector of length N containing the initial guess. (Input)
Default: XGUESS = 0.0.

RNORM — The Euclidean length of components of the function f (x) after the approximate solution has been found. (Output).

ISTAT — Scalar indicating further information about the approximate solution X. (Output)
See the Comments section for a description of the tolerances and the vectors IPARAM and RPARAM.

ISTAT Meaning

1 The function f (x) has a length less than TOLF = RPARAM(1). This is the expected value for ISTAT when an actual zero value of f (x) is anticipated.

2 The function f (x) has reached a local minimum. This is the expected value for ISTAT when a nonzero value of f (x) is anticipated.

3 A small change (absolute) was noted for the vector x. A full model problem step was taken. The condition for ISTAT = 2 may also be satisfied, so that a minimum has been found. However, this test is made before the test for ISTAT = 2.

4 A small change (relative) was noted for the vector x. A full model problem step was taken. The condition for ISTAT = 2 may also be satisfied, so that a minimum has been found. However, this test is made before the test for ISTAT = 2.

5 The number of terms in the quadratic model is being restricted by the amount of storage allowed for that purpose. It is suggested, but not required, that additional storage be given for the quadratic model parameters. This is accessed through the vector
IPARAM, documented below.

6 Return for evaluation of function and Jacobian if reverse
communication is desired. See the Comments below.

FORTRAN 90 Interface

Generic: CALL BCNLS (FCN, M, C, BL, BU, IRTYPE, XLB, XUB, X [,…])

Specific: The specific interface names are S_BCNLS and D_BCNLS.

FORTRAN 77 Interface

Single: CALL BCNLS (FCN, M, N, MCON, C, LDC, BL, BU, IRTYPE, XLB, XUB, XGUESS, X, RNORM, ISTAT)

Double: The double precision name is DBCNLS.

Description

The routine BCNLS solves the nonlinear least squares problem

subject to

BCNLS is based on the routine DQED by R.J. Hanson and F.T. Krogh. The section of BCNLS that approximates, using finite differences, the Jacobian of f(x) is a modification of JACBF by D.E. Salane.

Comments

1. Workspace may be explicitly provided, if desired, by use of B2NLS/DB2NLS. The reference is:

CALL B2NLS (FCN, M, N, MCON, C, LDC, BL, BU, IRTYPE, XLB, XUB, XGUESS, X, RNORM,ISTAT, IPARAM, RPARAM, JAC, F, FJ, LDFJ, IWORK, LIWORK, WORK, LWORK)

The additional arguments are as follows:

IPARAM — Integer vector of length six used to change certain default attributes of BCNLS. (Input).
If the default parameters are desired for B2NLS, set IPARAM(1) to zero. Otherwise, if any nondefault parameters are desired for IPARAM or RPARAM, the following steps should be taken before calling B2NLS:

CALL B7NLS (IPARAM, RPARAM)
Set nondefault values for IPARAM and RPARAM.
If double precision is being used, DB7NLS should be called instead. Following is a list of parameters and the default values.

IPARAM(1) = Initialization flag.

IPARAM(2) = ITMAX, the maximum number of iterations allowed.

Default: 75

IPARAM(3) = a flag that suppresses the use of the quadratic model in the inner loop. If set to one, then the quadratic model is never used. Otherwise use the quadratic model where appropriate. This option decreases the amount of workspace as well as the computing overhead required. A user may wish to determine if the application really requires the use of the quadratic model.

Default: 0

IPARAM(4) = NTERMS, one more than the maximum number of terms used in the quadratic model.

Default: 5

IPARAM(5) = RCSTAT, a flag that determines whether forward or reverse communication is used. If set to zero, forward communication through functions FCN and JAC is used. If set to one, reverse communication is used, and the dummy routines B10LS/DB10LS and B11LS/DB11LS may be used in place of FCN and JAC, respectively. When BCNLS returns with ISTAT = 6, arrays F and FJ are filled with f(x) and the Jacobian of f(x), respectively. BCNLS is then called again.

Default: 0

IPARAM(6) = a flag that determines whether the analytic Jacobian, as supplied in JAC, is used, or if a finite difference approximation is computed. If set to zero, JAC is not accessed and finite differences are used. If set to one, JAC is used to compute the Jacobian.

Default: 0

RPARAM — Real vector of length 7 used to change certain default attributes of
BCNLS.   (Input)
For the description of RPARAM, we make the following definitions:
FC                                           current value of the length of f (x)
FB                                           best value of length of f (x)
FL                                           value of length of f (x) at the previous step
PV                                           predicted value of length of f (x), after the step is taken, using the approximating model ɛ machine epsilon = amach(4).
The conditions |FB - PV| ≤ TOLSNR*FB and |FC - PV| ≤ TOLP*FB and
|FC - FL| ≤ TOLSNR*FB together with taking a full model step, must be satisfied before the condition ISTAT = 2 is returned. (Decreasing any of the values for TOLF, TOLD, TOLX, TOLSNR, or TOLP will likely increase the number of iterations required for convergence.)

RPARAM(1) = TOLF, tolerance used for stopping when FC ≤ TOLF.

RPARAM(2) = TOLX, tolerance for stopping when change to x values has length less than or equal to TOLX*length of x values.

RPARAM(3) = TOLD, tolerance for stopping when change to x values has length less than or equal to TOLD.

RPARAM(4) = TOLSNR, tolerance used in stopping condition ISTAT = 2.

Default: 1.E-5

RPARAM(5) = TOLP, tolerance used in stopping condition ISTAT = 2.

Default: 1.E-5

RPARAM(6) = TOLUSE, tolerance used to avoid values of x in the quadratic model's interpolation of previous points. Decreasing this value may result in more terms being included in the quadratic model.

RPARAM(7) = COND, largest condition number to allow when solving for the quadratic model coefficients. Increasing this value may result in more terms being included in the quadratic model.

Default: 30

JAC — User-supplied subroutine to evaluate the Jacobian. The usage is
CALL JAC(M, N, X, FJAC, LDFJAC), where
M - Number of functions.   (Input)
N - Number of variables.   (Input)
X - Array of length N containing the point at which the Jacobian will be evaluated.   (Input)
FJAC - The computed M ´ N Jacobian at the point X.   (Output)
LDFJAC - Leading dimension of the array FJAC.   (Input)
The routine JAC must be declared EXTERNAL in the calling program.

F — Real vector of length N used to pass f(x) if reverse communication
(IPARAM(4)) is enabled. (Input)

FJ — Real array of size M ´ N used to store the Jacobian matrix of f(x) if reverse communication (IPARAM(4)) is enabled. (Input)
Specifically,

LDFJ — Leading dimension of FJ exactly as specified in the dimension statement of the calling program. (Input)

IWORK — Integer work vector of length LIWORK.

LIWORK — Length of work vector IWORK. LIWORK must be at least
5MCON + 12N + 47 + MAX(M, N)

WORK — Real work vector of length LWORK

LWORK — Length of work vector WORK. LWORK must be at least
41N + 6M + 11MCON + (M + MCON)(N + 1) + NA(NA + 7) + 8 MAX(M, N) + 99.
Where NA = MCON + 2N + 6.

2. Informational errors

Type Code

3 1 The function f (x) has reached a value that may be a local minimum. However, the bounds on the trust region defining the size of the step are being hit at each step. Thus, the situation is suspect. (Situations of this type can occur when the solution is at infinity at some of the components of the unknowns, x).

3 2 The model problem solver has noted a value for the linear or quadratic model problem residual vector length that is greater than or equal to the current value of the function, i.e. the Euclidean length of f (x). This situation probably means that the evaluation of f (x) has more uncertainty or noise than is possible to account for in the tolerances used to not a local minimum. The value of x is suspect, but a minimum has probably been found.

3 3 More than ITMAX iterations were taken to obtain the solution. The value obtained for x is suspect, although it is the best set of x values that occurred in the entire computation. The value of ITMAX can be increased though the IPARAM vector.

Example 1

This example finds the four variables x1, x2, x3, x4 that are in the model function

There are values of h(t) at five values of t.

h(0.05) = 2.206

h(0.1) = 1.994

h(0.4) = 1.35

h(0.5) = 1.216

h(1.0) = 0.7358

There are also the constraints that x2, x4 ≤ 0, x1, x3 ≥ 0, and x2 and x4 must be separated by at least 0.05. Nothing more about the values of the parameters is known so the initial guess is 0.

USE BCNLS_INT

USE UMACH_INT

USE WRRRN_INT

IMPLICIT NONE

INTEGER MCON, N

PARAMETER (MCON=1, N=4)

! SPECIFICATIONS FOR PARAMETERS

INTEGER LDC, M

PARAMETER (M=5, LDC=MCON)

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER IRTYPE(MCON), NOUT

REAL BL(MCON), C(MCON,N), RNORM, X(N), XLB(N), &

XUB(N)

! SPECIFICATIONS FOR SUBROUTINES

! SPECIFICATIONS FOR FUNCTIONS

EXTERNAL FCN

CALL UMACH (2, NOUT)

! Define the separation between x(2)

! and x(4)

C(1,1) = 0.0

C(1,2) = 1.0

C(1,3) = 0.0

C(1,4) = -1.0

BL(1) = 0.05

IRTYPE(1) = 2

! Set lower bounds on variables

XLB(1) = 0.0

XLB(2) = 1.0E30

XLB(3) = 0.0

XLB(4) = 1.0E30

! Set upper bounds on variables

XUB(1) = -1.0E30

XUB(2) = 0.0

XUB(3) = -1.0E30

XUB(4) = 0.0

CALL BCNLS (FCN, M, C, BL, BL, IRTYPE, XLB, XUB, X, RNORM=RNORM)

CALL WRRRN ('X', X, 1, N, 1)

WRITE (NOUT,99999) RNORM

99999 FORMAT (/, 'rnorm = ', E10.5)

END

SUBROUTINE FCN (M, N, X, F)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER M, N

REAL X(*), F(*)

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER I

! SPECIFICATIONS FOR SAVE VARIABLES

REAL H(5), T(5)

SAVE H, T

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC EXP

REAL EXP

DATA T/0.05, 0.1, 0.4, 0.5, 1.0/

DATA H/2.206, 1.994, 1.35, 1.216, 0.7358/

DO 10 I=1, M

F(I) = X(1)*EXP(X(2)*T(I)) + X(3)*EXP(X(4)*T(I)) - H(I)

10 CONTINUE

RETURN

END

Output

                   X
       1       2       3       4
   1.999 -1.000   0.500 -9.954

rnorm = .42425E-03

Additional Examples

Example 2

This example solves the same problem as the last example, but reverse communication is used to evaluate f(x) and the Jacobian of f(x). The use of the quadratic model is turned off.

USE B2NLS_INT

USE UMACH_INT

USE WRRRN_INT

IMPLICIT NONE

INTEGER LDC, LDFJ, M, MCON, N

PARAMETER (M=5, MCON=1, N=4, LDC=MCON, LDFJ=M)

! Specifications for local variables

INTEGER I, IPARAM(6), IRTYPE(MCON), ISTAT, IWORK(1000), &

LIWORK, LWORK, NOUT

REAL BL(MCON), C(MCON,N), F(M), FJ(M,N), RNORM, RPARAM(7), &

WORK(1000), X(N), XGUESS(N), XLB(N), XUB(N)

REAL H(5), T(5)

SAVE H, T

INTRINSIC EXP

REAL EXP

! Specifications for subroutines

EXTERNAL B7NLS

! Specifications for functions

EXTERNAL B10LS, B11LS

DATA T/0.05, 0.1, 0.4, 0.5, 1.0/

DATA H/2.206, 1.994, 1.35, 1.216, 0.7358/

CALL UMACH (2, NOUT)

! Define the separation between x(2)

! and x(4)

C(1,1) = 0.0

C(1,2) = 1.0

C(1,3) = 0.0

C(1,4) = -1.0

BL(1) = 0.05

IRTYPE(1) = 2

! Set lower bounds on variables

XLB(1) = 0.0

XLB(2) = 1.0E30

XLB(3) = 0.0

XLB(4) = 1.0E30

! Set upper bounds on variables

XUB(1) = -1.0E30

XUB(2) = 0.0

XUB(3) = -1.0E30

XUB(4) = 0.0

! Set initial guess to 0.0

XGUESS = 0.0E0

! Call B7NLS to set default parameters

CALL B7NLS (IPARAM, RPARAM)

! Suppress the use of the quadratic

! model, evaluate functions and

! Jacobian by reverse communication

IPARAM(3) = 1

IPARAM(5) = 1

IPARAM(6) = 1

LWORK = 1000

LIWORK = 1000

! Specify dummy routines for FCN

! and JAC since we are using reverse

! communication

10 CONTINUE

CALL B2NLS (B10LS, M, N, MCON, C, LDC, BL, BL, IRTYPE, XLB, &

XUB, XGUESS, X, RNORM, ISTAT, IPARAM, RPARAM, &

B11LS, F, FJ, LDFJ, IWORK, LIWORK, WORK, LWORK)

! Evaluate functions if the routine

! returns with ISTAT = 6

IF (ISTAT .EQ. 6) THEN

DO 20 I=1, M

FJ(I,1) = EXP(X(2)*T(I))

FJ(I,2) = T(I)*X(1)*FJ(I,1)

FJ(I,3) = EXP(X(4)*T(I))

FJ(I,4) = T(I)*X(3)*FJ(I,3)

F(I) = X(1)*FJ(I,1) + X(3)*FJ(I,3) - H(I)

20 CONTINUE

GO TO 10

END IF

CALL WRRRN ('X', X, 1, N, 1)

WRITE (NOUT,99999) RNORM

99999 FORMAT (/, 'rnorm = ', E10.5)

END

Output

                   X
       1       2       3       4
   1.999 -1.000   0.500 -9.954

rnorm = .42413E-03

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260