/* Program to accompany the article "How to find a feasible point for a constrained optimization in SAS" by Rick Wicklin, published 19JUN2017 on The DO Loop blog, http://blogs.sas.com/content/iml/2017/06/19/find-feasible-point-optimization-sas.html */ proc iml; con = { 0 0 . ., /* param min */ 10 10 . ., /* param max */ 3 -2 -1 10, /* 3*x1 + -2*x2 LE 10 */ 5 10 -1 56, /* 5*x1 + 10*x2 LE 56 */ 4 2 1 7 }; /* 4*x1 + 2*x2 GE 7 */ guess = {0 0}; /* arbitrary p-dimensional point */ call nlpfea(z, guess, con); /* x0 is feasible point */ print guess[c={"x" "y"}], z[c={"Tx" "Ty"}]; start func(x); x1 = x[,1]; x2 = x[,2]; return ( -(x1-3)##4 + -(x2-2)##2 + 0.1*sin(x1#x2)); finish; opt = {1, /* find maximum of function */ 3}; /* print a little bit of output */ x0 = {0 0}; call nlpnra(rc, x_Opt, "func", x0, opt) blc=con; print rc; /* write feasible points to SAS data set */ NumPts = 36; twopi = 2*constant('pi'); x=.; y=.; Tx=.; Ty=.; create feasible var {"x" "y" "Tx" "Ty"}; do i = 1 to NumPts; x = 2.5 + 5*cos((i-1)/NumPts * twopi); /* guess on circle */ y = 2.5 + 5*sin((i-1)/NumPts * twopi); x = 3 + 4.2*cos((i-1)/NumPts * twopi); /* guess on circle */ y = 3 + 4.2*sin((i-1)/NumPts * twopi); call nlpfea(feasPt, x||y, con); /* transform into feasible */ Tx = feasPt[1]; Ty = feasPt[2]; append; end; close; QUIT; /* visualize the feasible region and feasible points */ data polygon; ID = 1; input PolyX PolyY; datalines; 1.75 0 3.3333 0 5.3 2.95 0 5.6 0 3.5 ; data All; set polygon feasible; run; ods graphics / width=500px height=480px; title "Transform Guess into the Feasible Region"; proc sgplot data=All aspect=1; polygon id=id x=PolyX y=PolyY /lineattrs=(thickness=3) fill outline name="F" legendlabel="Feasible Region"; scatter x=x y=y / markerattrs=(symbol=CircleFilled size=6) name="G" legendlabel="Guess"; scatter x=Tx y=Ty / markerattrs=GraphData2(symbol=CircleFilled size=10) name="P" legendlabel="Transformed" ; vector x=Tx y=Ty / xorigin=x yorigin=y lineattrs=(color=black); xaxis display=(nolabel); yaxis display=(nolabel); keylegend "F" "G" "P"; run;