Quadratic Programming

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: QUADPROG                                            */
/*   TITLE: Quadratic Programming                               */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS: MATRIX  OR      SUGI6                               */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: MDC                         UPDATE:                 */
/*     REF:                                                     */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/


proc iml;
start qp( names, c, H, G, rel, b, activity);
   if min(eigval(h))<0 then do;
      error={'The minimum eigenvalue of the H matrix is negative.',
             'Thus it is not positive semidefinite.',
             'QP is terminating.'};
      print error;
      stop;
   end;
   nr=nrow(G);
   nc=ncol(G);

   /* Put in canonical form */
   rev = (rel='<=');
   adj = (-1 * rev) + ^rev;
   g = adj# G;
   b = adj # b;
   eq = ( rel = '='  );
   if max(eq)=1 then do;
      g = g // -(diag(eq)*G)[loc(eq),];
      b = b // -(diag(eq)*b)[loc(eq)];
   end;
   m = (h || -g`) // (g || j(nrow(g),nrow(g),0));
   q = c // -b;

   /* Solve the problem */
   call lcp(rc,w,z,M,q);

   /* Report the solution */
   print ({'*************Solution is optimal***************',
           '*********No solution possible******************',
           ' ', ' ', ' ',
           '**********Solution is numerically unstable*****',
           '***********Not enough memory*******************',
           '**********Number of iterations exceeded********'}[rc+1]);
   activity = z[1:nc];
   objval = c`*activity + activity`*H*activity/2;
   print objval[L='Objective Value'],
         activity[r=names L= 'Decision Variables'];
finish qp;

c = { 0, 0, 0, 0 };
h = { 1003.1 4.3 6.3 5.9 ,
         4.3 2.2 2.1 3.9 ,
         6.3 2.1 3.5 4.8 ,
         5.9 3.9 4.8 10 };
g = {  1   1   1   1  ,
      .17 .11 .10 .18 };

/* constraints: proportions sum to 1; gain at least 10% */
b = { 1 , .10 };
rel = { '=', '>='};
names = 'Asset1':'Asset4';
run qp(names, c, h, g, rel, b, activity);