``` /******************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                   */
/*                                                                */
/*   TITLE: Quadratic Programming Algorithm                       */
/* PRODUCT: IML                                                   */
/*  SYSTEM: ALL                                                   */
/*    KEYS: MATH MATRIX                                           */
/*   PROCS: IML                                                   */
/*    DATA:                                                       */
/*                                                                */
/* SUPPORT: RHD                         UPDATE:                   */
/*     REF: METHOD OF THEIL & VAN DE PANNE                        */
/*    MISC: CONVERTED FROM MATRIX TO IML USING MATIML             */
/******************************************************************/

*---------------- QUADRATIC PROGRAMMING IN IML ------------------*
| THE FOLLOWING MATRIX PROGRAM IMPLEMENTS THE METHOD OF THEIL &  |
| VAN DE PANNE (MANAGEMENT SCIENCE V.7,PP1-20 (1960)) FOR        |
| MAXIMIZING A QUADRATIC FUNCTION SUBJECT TO LINEAR CONSTRAINTS. |
| THE METHOD STARTS WITH THE OPTIMAL UNCONSTRINED VALUES AND THEN|
| SEARCHES VARIOUS CONFIGURATIONS OF THE CONSTRAINTS FOR THE     |
| OPTIMAL FEASIBLE VALUES.                                       |
|                                                                |
| THE OBVIOUS APPLICATION OF THIS METHOD IS IN LEAST SQUARES     |
| REGRESSION WITH INEQUALITY CONSTRAINTS ON THE PARAMETERS.      |
*---------------------------------------------------JPS 23OCT79--*;
PROC IML;
RESET AUTONAME LINESIZE=132 PAGESIZE=60;

START MAIN;

*---DATA FOR THIS PROBLEM---;
XPY={18,16,22,20};
XPX={ 6  1  8  0,
1 10  1  4,
8  1 17  3,
0  4  3 11};
R=  {-1  0  0  0,
0 -1  0  0,
0  0 -1  0,
0  0  0 -1,
1  1  1  1,
5  0 10  0,
0  4  0  5};
C={0 0 0 0 1.66666667 2 3}`;

*---FIND DIMENSIONS OF PROBLEM---;
NX= NROW(XPX);                  * NUMBER OF VARIABLES;
NR= NROW(R);                    * NUMBER OF INEQ. CONSTRAINTS;
NT=NX+NR+1;                     * TOTAL COLUMNS IN Z;
IR=NX+1:NX+NR;                  * INDEX FOR CONSTRAINTS;
Z1= J(NR,NR,0);

*---SET UP CROSSPRODUCTS---;
Z= ( XPX || R` || XPY )//
(  R  || Z1 ||  C  )//
(XPY` || C` ||  0  );

*---SOLVE FOR UNCONSTRAINED MAXIMUM---;
Z= SWEEP(Z,1:NX);
ZTRY=Z;
I=1;
NS=1;

*---LOOK FOR UNSATISFIED CONSTRAINTS---;
STABLE= Z[IR,NT]<0;
CTABLE= J(NR,1,0);
IF ( SUM(STABLE)=0) THEN GOTO DONE;

LOOP: C= CTABLE[,I];
S= STABLE[,I];
NLOC= LOC(S);
NNEG= NCOL(NLOC);
DO J=1 TO NNEG;
IJ= NLOC[,J];
CTRY=C;
CTRY[IJ,]=1;
* FORM NEW COMBINATION;
DO K=I+1 TO NS;
* SEE IF ALREADY IN STABLE;
IF ( CTRY= CTABLE[,K]) THEN GOTO DUP;
END;
CTABLE=CTABLE||CTRY;
NS=NS+1;

STABLE=STABLE||STRY;
IF ( SUM(STRY)=0) THEN  LINK TEST;

DUP: END;

I=I+1;
IF ( I<=NS) THEN  GOTO LOOP;
PRINT " -----NO FEASIBLE SOLUTION-----";
STOP;

*---ROUTINE TO SWEEP IN CONSTRAINTS---;
SWP: ITRY= LOC(CTRY)+NX;
ZTRY= SWEEP(Z,ITRY);
STRY= ZTRY[IR,NT]<0;
STRY=STRY#^CTRY;
RETURN;

*---ROUTINE TO TEST IF OPTIMAL---;
TEST: CTLOC= LOC(CTRY);
NTLOC= NCOL(CTLOC);
DO III=1 TO NTLOC;
IIJ= CTLOC[,III];
CT=CTRY;
CT[IIJ,]=0;
LOCCT= LOC(CT);
IF ( NCOL(LOCCT)>0) THEN  ZCT= SWEEP(Z,LOCCT+NX);
ELSE  ZCT=Z;
CJ=( ZCT[IR,NT] <0)#^CT;
IF ( ^CJ[IIJ,]) THEN  RETURN;
END;

*---PRINT OUT RESULTS---;
DONE: PRINT CTRY[L="CONSTRAINT CONFIGURATION FOR FEASIBLE OPTIMUM"];
BETA= ZTRY[1:NX,NT];
PRINT BETA[L="PARAMETER ESTIMATES FOR OPTIMAL"];
LAMBDA= ZTRY[IR,NT];
PRINT LAMBDA[L="CONSTRAINT SLACKS AND LAGRANGIANS"];
PRINT CTABLE[F=4. L="CONFIGURATIONS TRIED"];
PRINT STABLE[F=4. L="PATTERN OF CONSTRAINT VIOLATIONS FOR CONFIGURATIONS"];
STOP;

FINISH MAIN;

RUN MAIN;

```