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;