Time-Optimal Heat Conduction

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: IMNLPEX8                                            */
/*   TITLE: Time-Optimal Heat Conduction                        */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswmh/sasghg               UPDATE: frwick 5/2015   */
/*     REF: Technical Report Examples                           */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

proc iml;
/* Vector mu[30] found by solving mu[j] * tan(mu[j]) = 1 */
mu = {
   8.6033358901938E-01 ,    3.4256184594817E+00 ,
   6.4372981791719E+00 ,    9.5293344053619E+00 ,
   1.2645287223856E+01 ,    1.5771284874815E+01 ,
   1.8902409956860E+01 ,    2.2036496727938E+01 ,
   2.5172446326646E+01 ,    2.8309642854452E+01 ,
   3.1447714637546E+01 ,    3.4586424215288E+01 ,
   3.7725612827776E+01 ,    4.0865170330488E+01 ,
   4.4005017920830E+01 ,    4.7145097736761E+01 ,
   5.0285366337773E+01 ,    5.3425790477394E+01 ,
   5.6566344279821E+01 ,    5.9707007305335E+01 ,
   6.2847763194454E+01 ,    6.5988598698490E+01 ,
   6.9129502973895E+01 ,    7.2270467060309E+01 ,
   7.5411483488848E+01 ,    7.8552545984243E+01 ,
   8.1693649235601E+01 ,    8.4834788718042E+01 ,
   8.7975960552493E+01 ,    9.1117161394464E+01  };

/* Vector a[nmu] depends only on mu and is computed once */
nmu  = nrow(mu);
a = j(1,nmu,0.);
do i = 1 to nmu;
   a[i] = 2*sin(mu[i]) / (mu[i] + sin(mu[i])*cos(mu[i]));
end;

/* This is the integrand used in h(x) */
start fquad(s) global(mu,rho);
   z = (rho * cos(s*mu) - 0.5*(1. - s##2))##2;
   return(z);
finish;

/* Obtain nonlinear constraint h(x) */
start h(x) global(n,nmu,mu,a,rho);
   xx = x##2;
   do i= n-1 to 1 by -1;
      xx[i] = xx[i+1] + xx[i];
   end;
   rho = j(1,nmu,0.);
   do i=1 to nmu;
      mu2 = mu[i]##2;
      sum = 0; t1n = -1.;
      do j=2 to n;
         t1n = -t1n;
         sum = sum + t1n * exp(-mu2*xx[j]);
      end;
      sum = -2*sum + exp(-mu2*xx[1]) + t1n;
      rho[i] = -a[i] * sum;
   end;
   aint = do(0,1,.5);
   call quad(z,"fquad",aint) eps=1.e-10;
   v = sum(z);
   return(v);
finish;

/* Define modules for NLPQN call: f, g, and c */
start F_HS88(x);
   f = x * x`;
   return(f);
finish F_HS88;

start G_HS88(x);
   g = 2 * x;
   return(g);
finish G_HS88;

start C_HS88(x);
   c = 1.e-4 - h(x);
   return(c);
finish C_HS88;

title 'Hock & Schittkowski Problem #91 (1981) n=5, INSTEP=1';
opt  = j(1,10,.);
opt[2]=3;
opt[10]=1;
tc   = j(1,12,.);
tc[6]=1.e-4;
x0 = {.5 .5 .5 .5 .5};
n = ncol(x0);

call nlpqn(rc,rx,"F_HS88",x0,opt,,tc) grd="G_HS88" nlc="C_HS88";
quit;