Chemical Equilibrium

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: IMNLPEX1                                            */
/*   TITLE: Chemical Equilibrium                                */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswmh                      UPDATE: frwick 5/2015   */
/*     REF: Technical Report Examples                           */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

proc iml;
c = { -6.089 -17.164 -34.054  -5.914 -24.721
      -14.986 -24.100 -10.708 -26.662 -22.179 };
start F_BRACK(x) global(c);
   s = x[+];
   f = sum(x # (c + log(x / s)));
   return(f);
finish F_BRACK;

con = { .  .  .  .  .  .  .  .  .  .    .   .  ,
        .  .  .  .  .  .  .  .  .  .    .   .  ,
        1. 2. 2. .  .  1. .  .  .  1.   0.  2. ,
        .  .  .  1. 2. 1. 1. .  .  .    0.  1. ,
        .  .  1. .  .  .  1. 1. 2. 1.   0.  1. };
con[1,1:10] = 1.e-6;

x0 = j(1,10, .1);
optn = {0 3};

title 'NLPTR subroutine: No Derivatives';
call nlptr(xres,rc,"F_BRACK",x0,optn,con);

options nodate ls=68 ps=60;
proc iml;
c = { -6.089 -17.164 -34.054  -5.914 -24.721
    -14.986 -24.100 -10.708 -26.662 -22.179 };

start F_BRACK(x) global(c);
   s = x[+];
   f = sum(x # (c + log(x / s)));
   return(f);
finish F_BRACK;

start G_BRACK(x) global(c);
   s = x[+];
   g = c + log(x / s);
   return(g);
finish G_BRACK;

con = { .  .  .  .  .  .  .  .  .  .    .   .  ,
      .  .  .  .  .  .  .  .  .  .    .   .  ,
      1. 2. 2. .  .  1. .  .  .  1.   0.  2. ,
      .  .  .  1. 2. 1. 1. .  .  .    0.  1. ,
      .  .  1. .  .  .  1. 1. 2. 1.   0.  1. };
con[1,1:10] = 1.e-6;

x0 = j(1,10, .1);
optn = {0 3};

CALL NLPTR(xres,rc,"F_BRACK",x0,optn,con);