*---------------------------------------------------------
* Date created:  February 25, 2001
* 
* Brown, A.A., Bartholomew-Biggs,M.C.,  
* ODE versus SQP methods for constrained optimisation
* Technical Report No. 179, June, 1987, 
* The Hatfield Polytechnic, College Lane, Hatfield,
* Problem 14.     
*
* Trajectory optimization problem     
*
* Dr. Neculai Andrei,
* Research Institute for Informatics - Bucharest
* 8-10, Averescu Avenue,
* 71316 Bucharest - Romania      
* E-mail: nandrei@u3.ici.ro   
* web: http://www.ici.ro/camo/neculai/nandrei.htm
*---------------------------------------------------------
*
      subroutine ini(n,m,mb,me,sb,x,icgc,ipgc,icge,ipge,
     1		     nszc,nsze,nb)
*
      double precision x(n)
      integer sb(mb)
      integer icgc(nszc),ipgc(n+1)
      integer icge(nsze),ipge(n+1)     
*
	write(1,10)
10	format(5x,'Example of a Nonlinear Programming Problem.')
	write(1,11)
11	format(5x,'Brown - Bartholomew-Biggs, 1987, Problem 14.')  
      write(1,12)
12    format(5x,'Trajectory optimization problem')      
*
      n=8
      m=0
      me=6
      mb=16
      nszc=0
      nsze=21
*
* Initial point
*
      x(1)=5.67d0
      x(2)=5.23d0
      x(3)=11.96d0
      x(4)=23.88d0
      x(5)=-2.61d0
      x(6)=2.34d0
      x(7)=3.45d0
      x(8)=4.89d0
*
* Index vector for simple bounds:
*              
      j=1
      do i=1,n
	  sb(j)=i
	  sb(j+1)=-i
	  j=j+2
      end do
*
* Rows indices of the nonzeros of the Jacobian of the Equalities.
*               
      icge(1)=1
      icge(2)=4
      icge(3)=5
      icge(4)=6
      
      icge(5)=4
      icge(6)=5
      icge(7)=6
      
      icge(8)=2
      icge(9)=4 
      icge(10)=5
      icge(11)=6
      
      icge(12)=3
      icge(13)=4
      icge(14)=5
      icge(15)=6
      
      icge(16)=4
      icge(17)=5
      icge(18)=6
      
      icge(19)=1
      icge(20)=2
      icge(21)=3
*
* Starting address of the columns of the Jacobian of the Equalities.
*
      ipge(1)=1
      ipge(2)=5
      ipge(3)=8
      ipge(4)=12
      ipge(5)=16  
      ipge(6)=19
      ipge(7)=20
      ipge(8)=21
      ipge(9)=22
*
* Rows indices of the nonzeros of the Jacobian of the INEqualities.
*                                                          
*
* Starting address of the columns of the Jacobian of the INEqualities.
*             
*
      return
      end
*


*--------------------------------------------------------------
* Date created:  February 25, 2001
* Brown - Bartholomew-Biggs, 1987, Problem 14.         
* Trajectory optimization problem
*--------------------------------------------------------------
*
      subroutine prob(n,m,mb,me,sb,x,objf,gobj,c,gc,cb,e,ge,
     1                nszc,nsze,nb)
*
* Calculate problem functions at iterate x.
*
      double precision x(n),objf,gobj(n),c(m),gc(nszc),
     *   cb(mb),e(me),ge(nsze)
*
      real*8 z,y,u,v
      real*8 dzx1, dzx2, dzx3, dzx4, dzx5
      real*8 dyx1, dyx2, dyx3, dyx4, dyx5
      real*8 dux1, dux2, dux4, dux5
      real*8 dvx1, dvx2, dvx4, dvx5     
*
* Objective function, and its gradient.
*
      objf=x(1) + x(3)/2.d0 + x(4)
*
	gobj(1)=1.d0
	gobj(2)=0.d0
	gobj(3)=0.5d0
	gobj(4)=1.d0
	gobj(5)=0.d0
	gobj(6)=0.d0
	gobj(7)=0.d0
	gobj(8)=0.d0
*
* Bounds on variables.
*
      j=1
	do i=1,n
	  cb(j)=x(i) + 1.d38  
	  cb(j+1)=1.d38 - x(i)
	  j=j+2
	end do       
*
* INEquality Constraints.
*
*
* Jacobian of the INEquality constraints.
*
*
* Equality constraints
*          
      z=0.1d0*(x(1)+x(3)+x(4)) + 0.01d0*dcos(x(2))*
     *        (x(1)**2+2.d0*x(1)*x(3)+2.d0*x(1)*x(4)) +
     *        0.01d0*x(4)*x(4)*dcos(x(5))
*     
      y=1.d0 + 0.1d0*(x(1)+x(3)+x(4)) + 0.01d0*dsin(x(2))*
     *        (x(1)**2+2.d0*x(1)*x(3)+2.d0*x(1)*x(4)) +
     *        0.01d0*x(4)*x(4)*dsin(x(5))
*     
      u=0.1d0 + 0.02d0*( x(1)*dcos(x(2))+x(4)*dcos(x(5)) )
      v=0.1d0 + 0.02d0*( x(1)*dsin(x(2))+x(4)*dsin(x(5)) )
*               
*-----------
*
      e(1)=x(1) - x(6)*x(6) - 0.01d0
*      
      e(2)=x(3) - x(7)*x(7)
*      
      e(3)=x(4) - x(8)*x(8)
*      
      e(4)=z*z + y*y - 4.0d0
*      
      e(5)=z*u + v*y
*      
      e(6)=u*y - v*z - 0.7d0 
*      
* Jacobian of the Equality constraints
*
*dz 
      dzx1=0.1d0 + 0.01d0*dcos(x(2))*( 2.d0*x(1)+2.d0*x(3)+2.d0*x(4) )
      dzx2=0.01d0*(x(1)**2+2.d0*x(1)*x(3)+2.d0*x(1)*x(4))*(-dsin(x(2)))
      dzx3=0.1d0+0.01d0*dcos(x(2))*2.d0*x(1)
      dzx4=0.1d0+0.01d0*dcos(x(2))*2.d0*x(1) + 0.02d0*x(4)*dcos(x(5))
      dzx5=-0.01d0*x(4)*x(4)*dsin(x(5))
*dy
      dyx1=0.1d0 + 0.01d0*dsin(x(2))*( 2.d0*x(1)+2.d0*x(3)+2.d0*x(4) )
      dyx2=0.01d0*(x(1)**2+2.d0*x(1)*x(3)+2.d0*x(1)*x(4))*dcos(x(2))
      dyx3=0.1d0+0.01d0*dsin(x(2))*2.d0*x(1)
      dyx4=0.1d0+0.01d0*dsin(x(2))*2.d0*x(1) + 0.02d0*x(4)*dsin(x(5))
      dyx5=0.01d0*x(4)*x(4)*dcos(x(5))
*du
      dux1= 0.02d0*dcos(x(2))
      dux2=-0.02d0*x(1)*dsin(x(2))      
      dux4= 0.02d0*dcos(x(5))
      dux5=-0.02d0*x(4)*dsin(x(5))
*dv
      dvx1=0.02d0*dsin(x(2))
      dvx2=0.02d0*x(1)*dcos(x(2))
      dvx4=0.02d0*dsin(x(5))
      dvx5=0.02d0*x(4)*dcos(x(5))      
*      
c1                        
      ge(1)=1.d0
      ge(2)=2.d0*z*dzx1 + 2.d0*y*dyx1
      ge(3)=z*dux1 + u*dzx1 + v*dyx1 + y*dvx1
      ge(4)=u*dyx1 + y*dux1 - v*dzx1 - z*dvx1
c2     
      ge(5)=2.d0*z*dzx2 + 2.d0*y*dyx2
      ge(6)=z*dux2 + u*dzx2 + v*dyx2 + y*dvx2
      ge(7)=u*dyx2 + y*dux2 - v*dzx2 - z*dvx2
c3     
      ge(8)= 1.d0
      ge(9)= 2.d0*z*dzx3 + 2.d0*y*dyx3
      ge(10)=z*dux3 + u*dzx3 + v*dyx3 + y*dvx3
      ge(11)=u*dyx3 + y*dux3 - v*dzx3 - z*dvx3
c4      
      ge(12)=1.d0
      ge(13)=2.d0*z*dzx4 + 2.d0*y*dyx4
      ge(14)=z*dux4 + u*dzx4 + v*dyx4 + y*dvx4 
      ge(15)=u*dyx4 + y*dux4 - v*dzx4 - z*dvx4
c5     
      ge(16)=2.d0*z*dzx5 + 2.d0*y*dyx5
      ge(17)=z*dux5 + u*dzx5 + v*dyx5 + y*dvx5 
      ge(18)=u*dyx5 + y*dux5 - v*dzx5 - z*dvx5 
c6     
      ge(19)=-2.d0*x(6)
c7      
      ge(20)=-2.d0*x(7)
c8
      ge(21)=-2.d0*x(8)      
*      
      return
      end
*---------------------------------------------bb14.for