Try out NAG Library functions

Explore NAG maths and stats routines with interactive demos
Function ID
E05USF
Name
nagf_glopt_nlp_multistart_sqp_lsq
Description
Global optimization of a sum of squares problem using multi-start, nonlinear constraints
Keywords
multi-start algorithm | sequential QP method | sum of squares
This example is based on Problem 57 in Hock and Schittkowski (1981) and involves the minimization of the sum of squares function
Fx=12i=144yi-fix2,  
where
fix=x1+0.49-x1e-x2ai-8  
and
iyiaiiyiai10.498230.412220.498240.402230.4810250.422440.4710260.402450.4810270.402460.4710280.412670.4612290.402680.4612300.412690.4512310.4128100.4312320.4028110.4514330.4030120.4314340.4030130.4314350.3830140.4416360.4132150.4316370.4032160.4316380.4034170.4618390.4136180.4518400.3836190.4220410.4038200.4220420.4038210.4320430.3940220.4122440.3942  
subject to the bounds
x1-0.4x2-4.0  
to the general linear constraint
x1+x21.0  
and to the nonlinear constraint
0.49x2-x1x20.09.  
The example data reflects that shown in the "Example" section of the routine documentation. You can change this here to try alternative inputs. The formatting will need to be kept as it is here, otherwise the program is likely to fail to run correctly.

Please note that incompatible data will however cause the example output to display an error message. These error messages are fully explained in the Routine document
!   E05USF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module e05usfe_mod

!     E05USF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: confun, objfun, start
    Contains
      Subroutine objfun(mode,m,n,ldfjsl,needfi,x,f,fjsl,nstate,iuser,ruser)
!       Evaluates the subfunctions and their 1st derivatives.

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: a(44) = (/8._nag_wp,8._nag_wp,       &
                                          10._nag_wp,10._nag_wp,10._nag_wp,    &
                                          10._nag_wp,12._nag_wp,12._nag_wp,    &
                                          12._nag_wp,12._nag_wp,14._nag_wp,    &
                                          14._nag_wp,14._nag_wp,16._nag_wp,    &
                                          16._nag_wp,16._nag_wp,18._nag_wp,    &
                                          18._nag_wp,20._nag_wp,20._nag_wp,    &
                                          20._nag_wp,22._nag_wp,22._nag_wp,    &
                                          22._nag_wp,24._nag_wp,24._nag_wp,    &
                                          24._nag_wp,26._nag_wp,26._nag_wp,    &
                                          26._nag_wp,28._nag_wp,28._nag_wp,    &
                                          30._nag_wp,30._nag_wp,30._nag_wp,    &
                                          32._nag_wp,32._nag_wp,34._nag_wp,    &
                                          36._nag_wp,36._nag_wp,38._nag_wp,    &
                                          38._nag_wp,40._nag_wp,42._nag_wp/)
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ldfjsl, m, n, needfi, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(m)
        Real (Kind=nag_wp), Intent (Inout) :: fjsl(ldfjsl,*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: temp
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..

!       This is a two-dimensional objective function.
!       As an example of using the mode mechanism,
!       terminate if any other problem size is supplied.

        If (n/=2) Then
          mode = -1
        End If

        If (nstate==1) Then
!         This is the first call.
!         Take any special action here if desired.
          Continue
        End If

        If (mode==0 .And. needfi>0) Then
          f(needfi) = x(1) + (0.49_nag_wp-x(1))*exp(-x(2)*(a(needfi)-          &
            8.0_nag_wp))
        Else
          Do i = 1, m
            temp = exp(-x(2)*(a(i)-8._nag_wp))

            If (mode==0 .Or. mode==2) Then
              f(i) = x(1) + (0.49_nag_wp-x(1))*temp
            End If

            If (mode==1 .Or. mode==2) Then
              fjsl(i,1) = 1._nag_wp - temp
              fjsl(i,2) = -(0.49_nag_wp-x(1))*(a(i)-8._nag_wp)*temp
            End If

          End Do
        End If

        Return
      End Subroutine objfun
      Subroutine confun(mode,ncnln,n,ldcjsl,needc,x,c,cjsl,nstate,iuser,ruser)
!       Evaluates the nonlinear constraints and their 1st derivatives.

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ldcjsl, n, ncnln, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: c(ncnln)
        Real (Kind=nag_wp), Intent (Inout) :: cjsl(ldcjsl,*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: needc(ncnln)
!       .. Executable Statements ..

!       This problem has only one constraint.
!       As an example of using the mode mechanism,
!       terminate if any other size is supplied.

        If (ncnln/=1) Then
          mode = -1
        End If

        If (nstate==1) Then

!         First call to CONFUN.  Set all Jacobian elements to zero.
!         Note that this will only work when 'Derivative Level = 3'
!         (the default; see Section 11.1 of the E04USA document).

          cjsl(1:ncnln,1:n) = 0._nag_wp
        End If

        If (needc(1)>0) Then

          If (mode==0 .Or. mode==2) Then
            c(1) = -0.09_nag_wp - x(1)*x(2) + 0.49_nag_wp*x(2)
          End If

          If (mode==1 .Or. mode==2) Then
            cjsl(1,1) = -x(2)
            cjsl(1,2) = -x(1) + 0.49_nag_wp
          End If

        End If

        Return
      End Subroutine confun
      Subroutine start(npts,quas,n,repeat1,bl,bu,iuser,ruser,mode)

!       Sets the initial points.
!       A typical user-defined start procedure.

!       .. Use Statements ..
        Use nag_library, Only: g05kgf, g05saf
!       .. Scalar Arguments ..
        Integer, Intent (Inout)        :: mode
        Integer, Intent (In)           :: n, npts
        Logical, Intent (In)           :: repeat1
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: bl(n), bu(n)
        Real (Kind=nag_wp), Intent (Inout) :: quas(n,npts), ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: genid, i, ifail, lstate, subid
!       .. Local Arrays ..
        Integer, Allocatable           :: state(:)
!       .. Executable Statements ..
!       quas is pre-assigned to zero.
        If (repeat1) Then
          quas(1,1) = 0.4_nag_wp
          quas(2,2) = 1._nag_wp
        Else
!         Generate a non-repeatable spread of points between bl and bu.
          genid = 2
          subid = 53
          lstate = -1
          Allocate (state(lstate))
          ifail = 0
          Call g05kgf(genid,subid,state,lstate,ifail)
          Deallocate (state)
          Allocate (state(lstate))
          ifail = 0
          Call g05kgf(genid,subid,state,lstate,ifail)
          Do i = 1, npts
            ifail = 0
            Call g05saf(n,state,quas(1,i),ifail)
            quas(1:n,i) = bl(1:n) + (bu(1:n)-bl(1:n))*quas(1:n,i)
          End Do
          Deallocate (state)
        End If
!       Set mode negative to terminate execution for any reason.
        mode = 0
        Return
      End Subroutine start
    End Module e05usfe_mod
    Program e05usfe

!     E05USF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: dgemv, e05usf, e05zkf, nag_wp
      Use e05usfe_mod, Only: confun, objfun, start
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: liopts = 740, lopts = 485, m = 44,   &
                                          n = 2, nb = 1, nclin = 1, ncnln = 1, &
                                          nin = 5, nout = 6, npts = 3
      Integer, Parameter               :: sdfjac = n
      Integer, Parameter               :: lda = nclin
      Integer, Parameter               :: ldc = ncnln
      Integer, Parameter               :: ldcjac = ncnln
      Integer, Parameter               :: ldclda = n + nclin + ncnln
      Integer, Parameter               :: ldfjac = m
      Integer, Parameter               :: ldx = n
      Integer, Parameter               :: listat = n + nclin + ncnln
      Logical, Parameter               :: repeat1 = .True.
!     .. Local Scalars ..
      Integer                          :: i, ifail, j, k, l, sda, sdcjac
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), bl(:), bu(:), c(:,:),        &
                                          cjac(:,:,:), clamda(:,:), f(:,:),    &
                                          fjac(:,:,:), work(:), x(:,:), y(:)
      Real (Kind=nag_wp)               :: objf(nb), opts(lopts), ruser(1)
      Integer                          :: info(nb), iopts(liopts), iter(nb),   &
                                          iuser(1)
      Integer, Allocatable             :: istate(:,:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'E05USF Example Program Results'
      Flush (nout)

!     Skip heading in data file.
      Read (nin,*)

      If (nclin>0) Then
        sda = n
      Else
        sda = 1
      End If

      If (ncnln>0) Then
        sdcjac = n
      Else
        sdcjac = 0
      End If

      Allocate (a(lda,sda),bl(n+nclin+ncnln),bu(n+nclin+ncnln),y(m),c(ldc,nb), &
        cjac(ldcjac,sdcjac,nb),f(m,nb),fjac(ldfjac,sdfjac,nb),                 &
        clamda(ldclda,nb),istate(listat,nb),x(ldx,nb),work(max(1,nclin)))

      If (nclin>0) Then
        Read (nin,*)(a(i,1:sda),i=1,nclin)
      End If

      Read (nin,*) y(1:m)
      Read (nin,*) bl(1:(n+nclin+ncnln))
      Read (nin,*) bu(1:(n+nclin+ncnln))

!     Initialize the solver.

      ifail = 0
      Call e05zkf('Initialize = E05USF',iopts,liopts,opts,lopts,ifail)

!     Solve the problem.

      ifail = -1
      Call e05usf(m,n,nclin,ncnln,a,lda,bl,bu,y,confun,objfun,npts,x,ldx,      &
        start,repeat1,nb,objf,f,fjac,ldfjac,sdfjac,iter,c,ldc,cjac,ldcjac,     &
        sdcjac,clamda,ldclda,istate,listat,iopts,opts,iuser,ruser,info,ifail)

      Select Case (ifail)
      Case (0)
        l = nb
      Case (8)
        l = info(nb)
        Write (nout,99999) iter(nb)
      Case Default
        Go To 100
      End Select

loop: Do i = 1, l
        Write (nout,99998) i
        Write (nout,99997) info(i)
        Write (nout,99996) 'Varbl'
        Do j = 1, n
          Write (nout,99995) 'V', j, istate(j,i), x(j,i), clamda(j,i)
        End Do
        If (nclin>0) Then
          Write (nout,99996) 'L Con'

!         Below is a call to the level 2 BLAS routine DGEMV.
!         This performs the matrix vector multiplication A*X
!         (linear constraint values) and puts the result in
!         the first NCLIN locations of WORK.

          Call dgemv('N',nclin,n,1.0_nag_wp,a,lda,x(1,i),1,0.0_nag_wp,work,1)

          Do k = n + 1, n + nclin
            j = k - n
            Write (nout,99995) 'L', j, istate(k,i), work(j), clamda(k,i)
          End Do
        End If
        If (ncnln>0) Then
          Write (nout,99996) 'N Con'
          Do k = n + nclin + 1, n + nclin + ncnln
            j = k - n - nclin
            Write (nout,99995) 'N', j, istate(k,i), c(j,i), clamda(k,i)
          End Do
        End If
        Write (nout,99994) objf(i)
        Write (nout,99993)
        Write (nout,99992)(clamda(k,i),k=1,n+nclin+ncnln)

        If (l==1) Then
          Exit loop
        End If

        Write (nout,*)

        Write (nout,*)                                                         &
          ' ------------------------------------------------------ '

      End Do loop

100   Continue

99999 Format (1X,I20,'starting points converged')
99998 Format (/,1X,'Solution number',I16)
99997 Format (/,1X,'Local minimization exited with code',I5)
99996 Format (/,1X,A,2X,'Istate',3X,'Value',9X,'Lagr Mult',/)
99995 Format (1X,A,2(1X,I3),4X,F12.4,2X,F12.4)
99994 Format (/,1X,'Final objective value = ',1X,F12.4)
99993 Format (/,1X,'QP multipliers')
99992 Format (1X,F12.4)
    End Program e05usfe
The NAG Library
The world’s largest collection of robust, documented, tested and maintained numerical algorithms.
Learn more here or contact us for purchasing information