*--------------------------------------------------------- * Date created: May 30, 2001 * * Mordechai Shacham, * Numerical solution of constrained nonlinear algebraic * equations. * * International Journal for Numerical Methods in Engineering, * Vol.23, 1986, pp.1455-1481. * * Problem 1, pp.1463. * Hiebert's 3rd Chemical Equilibrium 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,'Hiebert - 3rd chemical equilibrium problem ') * n=10 m=0 me=0 mb=10 nszc=0 nsze=0 * * Initial Point * x(1)=2.d0 x(2)=5.d0 x(3)=40.d0 x(4)=10.d0 x(5)=0.00001d0 x(6)=0.00001d0 x(7)=0.00001d0 x(8)=0.00001d0 x(9)=0.00001d0 x(10)=0.00001d0 * * Index vector for simple bounds: * do i=1,n sb(i)=i end do * * Rows indices of the nonzeros of the Jacobian of the Equalities. * * * Starting address of the columns of the Jacobian of the Equalities. * * * 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: May 30, 2001 * Hiebert's 3rd Chemical Equilibrium 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 r, s real*8 t(10), s1,s2,s3 real*8 a6,a7,a8,a9,a10 * * Objective function, and its gradient. * r=40 C s=0.d0 do i=1,10 s=s+x(i) end do C t(1)= x(1)+x(4)-3.d0 t(2)= 2.d0*x(1)+x(2)+x(4)+x(7)+x(8)+x(9)+2.d0*x(10)-r t(3)= 2.d0*x(2)+2.d0*x(5)+x(6)+x(7)-8.d0 t(4)= 2.d0*x(3)+x(5)-4.d0*r t(5)= x(1)*x(5)-0.193d0*x(2)*x(4) t(6)= x(6)*sqrt(x(2))-0.002597d0*sqrt(x(2)*x(4)*s) t(7)= x(7)*sqrt(x(4))-0.003448d0*sqrt(x(1)*x(4)*s) t(8)= x(8)*x(4)-1.799d0*x(2)*s/100000.d0 t(9)= x(9)*x(4)-2.155d0*x(1)*sqrt(x(3)*s)/10000.d0 t(10)= x(10)*x(4)*x(4)-3.846d0*x(4)*x(4)*s/100000.d0 C objf=0.d0 do i=1,10 objf = objf + t(i)*t(i) end do * * a6=0.002597d0 a7=0.003448d0 a8=1.799d0/100000.d0 a9=2.155d0/10000.d0 a10=3.846d0/100000.d0 * s1=sqrt(x(2)*x(4)*s) s2=sqrt(x(1)*x(4)*s) s3=sqrt(x(3)*s) * gobj(1)= 2.d0*t(1) + 4.d0*t(2) + 2.d0*t(5)*x(5) + * 2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1)) + * 2.d0*t(7)*(-a7*x(4)*(s+x(1))/(2.d0*s2)) + * 2.d0*t(9)*(-a9*(s3 + x(1)*x(3)/(2.d0*s3)))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(2)= 2.d0*t(2)+4.d0*t(3)+2.d0*t(5)*(-0.193d0*x(4))+ * 2.d0*t(6)*(x(6)/(2.d0*sqrt(x(2)))- * a6*x(4)*(x(2)+s)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*s)+ * 2.d0*t(9)*(-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(3)= 2.d0*t(4)+2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*x(2))+ * 2.d0*t(9)*(-a9*x(1)*(s+x(3))/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(4)= 2.d0*t(1)+2.d0*t(2)+2.d0*t(5)*(-0.193d0*x(2))+ * 2.d0*t(6)*(-a6*x(2)*(s+x(4))/(2.d0*s1))+ * 2.d0*t(7)*(x(7)/(2.d0*sqrt(x(4)))- * a7*x(1)*(x(4)+s)/(2.d0*s2))+ * 2.d0*t(8)*x(8)+ * 2.d0*t(9)*(x(9)-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(2.d0*x(4)*x(10)-a10*2.d0*x(4)*s) gobj(5)= 4.d0*t(3)+2.d0*t(4)+2.d0*t(5)*x(1)+ * 2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*x(2))+ * 2.d0*t(9)*(-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(6)= 2.d0*t(3)+ * 2.d0*t(6)*(sqrt(x(2))-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*x(2))+ * 2.d0*t(9)*(-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(7)= 2.d0*t(2)+2.d0*t(3)+ * 2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(sqrt(x(4))-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*x(2))+ * 2.d0*t(9)*(-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(8)= 2.d0*t(2)+2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(x(4)-a8*x(2))+ * 2.d0*t(9)*(-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(9)= 2.d0*t(2)+ * 2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*x(2))+ * 2.d0*t(9)*(x(4)-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*(-a10*x(4)*x(4)) gobj(10)= 4.d0*t(2)+ * 2.d0*t(6)*(-a6*x(2)*x(4)/(2.d0*s1))+ * 2.d0*t(7)*(-a7*x(1)*x(4)/(2.d0*s2))+ * 2.d0*t(8)*(-a8*x(2))+ * 2.d0*t(9)*(-a9*x(1)*x(3)/(2.d0*s3))+ * 2.d0*t(10)*x(4)*x(4)*(1.d0-a10) * * Bounds on variables. * do i=1,n cb(i)=x(i) end do * * Equality Constraints. * * * Jacobian of the equality constraints. * * * INEquality Constraints. * * * Jacobian of the INEquality constraints. * return end *----------------------------------------------------------CEQ3.for