call nlpqn within the loop

[Xing_]

Four functions: var_ete_pacrt , var_ete_dec_pacrt, funcre_pacrt, and mmd_pacrt_int are defined as follows: 

start var_ete_pacrt(numclust, clustsz, numrepperd, corrwtnperd, corrintperd, corrwtnpart) global(stddev, numperiod);

churnrate=numrepperd/clustsz;
corrclustperd=(clustsz*corrwtnperd*corrintperd+(1-churnrate)*(1-corrwtnperd)*corrwtnpart)/(1+(clustsz-1)*corrwtnperd);
f=4*stddev**2/(numperiod*numclust*clustsz)*(1+(clustsz-1)*corrwtnperd)*(1+(numperiod-1)*corrclustperd);

return (f);
finish var_ete_pacrt;

start var_ete_dec_pacrt(decclustsz) global(corrpar, stddev, numperiod, numrepperd, costperclust, costperpart, costperoutc, budget);
corrwtnperd = corrpar[1];
corrintperd = corrpar[2];
corrwtnpart = corrpar[3];

if decclustsz>0 then churnrate = numrepperd/decclustsz;

decnumclust=budget/(costperclust+costperpart*(1+(numperiod-1)*churnrate)*decclustsz+costperoutc*numperiod*decclustsz);
corrclustperd=(decclustsz*corrwtnperd*corrintperd+(1-corrwtnperd)*corrwtnpart*(1-churnrate))/(1+(decclustsz-1)*corrwtnperd);
f=4*stddev**2/(numperiod*decnumclust*decclustsz)*(1+(decclustsz-1)*corrwtnperd)*(1+(numperiod-1)*corrclustperd);

*print corrwtnperd corrintperd corrwtnpart corrclustperd decclustsz churnrate decnumclust f;
return (f);
finish var_ete_dec_pacrt;

start funcre_pacrt(x) global(corrpar, numclust, clustsz, minnumclust, maxnumclust, minclustsz, maxclustsz, numperiod, numrepperd, costperclust, costperpart, costperoutc, budget);
corrpar=j(3, 1);
corrpar[1:3] = x; * copy x -> corrpar;

print numclust clustsz;

corrwtnperd = x[1];
corrintperd = x[2];
corrwtnpart = x[3];

churncost=costperpart*(numperiod-1)*numrepperd;

bound1=(budget/maxnumclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);
bound2=(budget/minnumclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);

print corrpar, x, bound1 bound2;
con=j(3, 3, .);

con[1, 1]=max(2, bound1, minclustsz);
con[2, 1]=min(bound2, maxclustsz);

con[3, 1]=1;
con[3, 2]=1;
con[3, 3]=1;

par=j(10, 1, .);
par[9, 1]=4;

x0=10;
optn=j(1, 10, .);

optn[1]=0;
call nlpqn(rc, decclustsz, "var_ete_dec_pacrt", x0, optn) blc=con par=par;

re=999;
decnumclust=budget/(costperclust+(costperpart+costperoutc*numperiod)*decclustsz+churncost);
minvar=var_ete_dec_pacrt(decclustsz);

temp_var=var_ete_pacrt(decnumclust, decclustsz, numrepperd, corrwtnperd, corrintperd, corrwtnpart);
var=var_ete_pacrt(numclust, clustsz, numrepperd, corrwtnperd, corrintperd, corrwtnpart);
if var>0 & minvar>=0 then re = minvar/var;

return (re);
finish funcre_pacrt;

start mmd_pacrt_int(corrwtnperdmin, corrwtnperdmax, corrintperdmin, corrintperdmax, corrwtnpartmin, corrwtnpartmax)
global(costperclust, costperpart, costperoutc, budget, churncost,
numseq_lcrt, numperiod, numrepperd, clustszstart, minnumclust, maxnumclust, minclustsz, maxclustsz);

if mod(minnumclust, numseq_lcrt)=0 then numcluststart=minnumclust;
else numcluststart=(int(minnumclust/numseq_lcrt)+1)*numseq_lcrt;

reminmax=0;
do numclust=numcluststart to maxnumclust by numseq_lcrt;
clustsz=clustszstart;
actcost=numclust*(costperclust+(costperpart+costperoutc*numperiod)*clustsz+churncost);

clustsztemp=(budget/numclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);
clustszend=min(maxclustsz, int(clustsztemp));

do while (actcost<=budget & clustsz<=clustszend);
print numclust clustsz;
con=j(2, 3, .);

con[1, 1]=corrwtnperdmin;
con[2, 1]=corrwtnperdmax;
con[1, 2]=corrintperdmin;
con[2, 2]=corrintperdmax;
con[1, 3]=corrwtnpartmin;
con[2, 3]=corrwtnpartmax;

par=j(10, 1, .);
par[9, 1]=4;

x0=j(1, 3, .);
x0[1, 1]=(corrwtnperdmin+corrwtnperdmax)/2;
x0[1, 2]=(corrintperdmin+corrintperdmax)/2;
x0[1, 3]=(corrwtnpartmin+corrwtnpartmax)/2;

optn = {0};
call nlpqn(rc, xres, "funcre_pacrt", x0, optn) blc=con par=par;
remin=funcre_pacrt(xres);
print xres remin reminmax;

if remin>reminmax then
do;
mmdcorrwtnperd=xres[1];
mmdcorrintperd=xres[2];
mmdcorrwtnpart=xres[3];
mmdnumclust=numclust;
mmdclustsz=clustsz;
mmdre=remin;
mmdactcost=actcost;
reminmax=remin;
*print reflag numclust clustsz actcost remin revar retau reomega reminmax;
end;

clustsz=clustsz+1;
actcost=numclust*(costperclust+(costperpart+costperoutc*numperiod)*clustsz+churncost);
end; *enddowhile;
end; *end numclust;

mmdchurnrate=numrepperd/mmdclustsz;
mmd=corrwtnperdmin||corrwtnperdmax||corrintperdmin||corrintperdmax||corrwtnpartmin||corrwtnpartmax||mmdcorrwtnperd||mmdcorrintperd||mmdcorrwtnpart||numseq_lcrt||numperiod||numrepperd||mmdnumclust||mmdclustsz||mmdchurnrate||mmdre||mmdactcost;

return (mmd);
finish mmd_pacrt_int;

the main code:

trteff=1;
stddev=1;

corrwtnperdmin=0.1;
corrintperdmin=0.05;
corrwtnpartmin=0.2;
corrwtnperdmax=0.2;
corrintperdmax=0.1;
corrwtnpartmax=0.8;

costperclust=1000;
costperpart=300;
costperoutc=40;
budget=200000;

numrepperd=2;

numseq_lcrt=2;

numperiod=6;
minclustsz=2;
maxclustsz=100;
minnumclust=2;
maxnumclust=100;

clustszstart=max(minclustsz, numrepperd);
churncost=costperpart*(numperiod-1)*numrepperd;

mmdpacrt=j(1, 17);
mmdpacrt[1, 1:17]=mmd_pacrt_int(corrwtnperdmin, corrwtnperdmax, corrintperdmin, corrintperdmax, corrwtnpartmin, corrwtnpartmax);

Here is log file:

ERROR: Matrix numclust has not been set to a value.

statement : PRINT at line 2186 column 1
traceback : module FUNCRE_PACRT at line 2186 column 1
module MMD_PACRT_INT at line 2186 column 1

ERROR: Execution error as noted previously. (rc=1000)

operation : NLPQN at line 2186 column 1
operands : *LIT1130, x0, optn, con, , par

*LIT1130 1 row 1 col (character, size 12)

funcre_pacrt

x0 1 row 3 cols (numeric)

0.125 0.11 0.35

optn 1 row 1 col (numeric)

0

con 2 rows 3 cols (numeric)

0.05 0.02 0.2
0.2 0.2 0.5

par 10 rows 1 col (numeric)

.
.
.
.
.
.
.
.
4
.

statement : CALL at line 2186 column 1
traceback : module MMD_PACRT_INT at line 2186 column 1

statement : CALL at line 2186 column 1
traceback : module MMD_PACRT_INT at line 2186 column 1

ERROR: Matrix mmdpacrt has not been set to a value.

statement : PRINT at line 2186 column 1

[Rick_SAS]

The error message tells you that the variable 'numclust' is not defined when you call  FUNCRE_PACRT. The 'numclust' variable is defined in the GLOBAL clause, so it must be a global variable.

The FUNCRE_PACRT function is called by the NLPQN call in the MMD_PACRT_INT module from inside a loop that uses 'numclust' as a looping variable. The looping variable is local to the FUNCRE_PACRT function. It is not global.

If you want the 'numclust' variable to be accessible from outside the MMD_PACRT_INT module, add it to the GLOBAL clause for that function.

[Xing_]

Thank you for the response. I found another issue when finding the minimum of function var_ete_dec_pacrt. Below is my SAS syntax. 

proc iml;

start var_ete_dec_pacrt(decclustsz) global(corrpar, stddev, numperiod, numrepperd, costperclust, costperpart, costperoutc, budget);

corrwtnperd = corrpar[1];
corrintperd = corrpar[2];
corrwtnpart = corrpar[3];

if decclustsz>0 then churnrate = numrepperd/decclustsz;

decnumclust=budget/(costperclust+costperpart*(1+(numperiod-1)*churnrate)*decclustsz+costperoutc*numperiod*decclustsz);
corrclustperd=(decclustsz*corrwtnperd*corrintperd+(1-corrwtnperd)*corrwtnpart*(1-churnrate))/(1+(decclustsz-1)*corrwtnperd);

f=4*stddev**2/(numperiod*decnumclust*decclustsz)*(1+(decclustsz-1)*corrwtnperd)*(1+(numperiod-1)*corrclustperd);

return (f);
finish var_ete_dec_pacrt;

stddev=1;

costperclust=3000;
costperpart=200;
costperoutc=50;
budget=400000;
numrepperd=1;
numperiod=6;

minnumclust=20;
maxnumclust=40;
minclustsz=2;
maxclustsz=100;

corrwtnperd = 0.05;
corrintperd = 0.02;
corrwtnpart = 0.35;

corrpar=j(3, 1);
corrpar[1]=corrwtnperd;
corrpar[2]=corrintperd;
corrpar[3]=corrwtnpart;

churncost=costperpart*(numperiod-1)*numrepperd;
bound1=(budget/maxnumclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);
bound2=(budget/minnumclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);

con=j(2, 1, .);

con[1, 1]=max(numrepperd, bound1, minclustsz);
con[2, 1]=min(bound2, maxclustsz);

print corrpar con;

x0=10;
optn=j(1, 10, .);

optn[1]=0;
call nlpqn(rc, decclustsz, "var_ete_dec_pacrt", x0, optn) blc=con;
minvar=var_ete_dec_pacrt(decclustsz);

testclustsz=17;
testvar=var_ete_dec_pacrt(testclustsz);

print decclustsz minvar testclustsz testvar;
quit;

Function nlpqn did not provide the minimum as expected, since the return value is 0.042288, larger than the variance of 0.0042276 when I set cluster size as 17.

decclustsz minvar testclustsz testvar

16.788052

0.0042288

17

0.0042276

[Rick_SAS]

>Function nlpqn did not provide the minimum as expected, since the return value is 0.042288, larger than the variance of 0.0042276 when I set cluster size as 17.

Two comments:
1. The NLP functions assume that the underlying objective function is continuous. It sounds like you are trying to choose from a discrete set of clusters? If so, you don't need NLP for that. You can just perform a 1-D search. 

2. Your objective function is very flat near the optimal value, which is why the NLP routine is returning 16.788 instead of 17. The objective values differ in the 6th decimal place. You can try to tighten some of the convergence criteria, but essentially the result is telling you that it doesn't matter if you choose 16, 17, or 18 clusters. The variance is essentially the same in all cases. 

You can visualize the objective function to see how flat it is:

nc = T( 2:25 );
var = j(nrow(nc), 1, .);
do i = 1 to nrow(nc);
   var[i] = var_ete_dec_pacrt(nc[i]);
end;
/* find index of min */
idx = var[>:<];
ncMin = nc[idx];
varMin = var[idx];
print ncMin  varMin;

title "Objective function":
call scatter(nc, var) other="refline 0.0042252/axis=y;";

[Xing_]

cluster size can be decimal in the objective function. 

Is there a way to set up a tight termination criterion or any other solution? 

[Rick_SAS]

Yes, I said that. See the doc at SAS Help Center: Termination Criteria

[Xing_]

I did check the manual and tried several ways. None of them work. 

tc=j(13, 1, .);
*tc[6]=0.000000000001;
tc[7]=10;

[Rick_SAS]

It works fine for me:

tc=j(13, 1, .);
tc[6]=1e-8;
call nlpqn(rc, decclustsz, "var_ete_dec_pacrt", x0, optn) blc=con tc=tc;

[Xing_]

yes, it works. Other errors pop up.

start var_ete_pacrt(numclust, clustsz, numrepperd, corrwtnperd, corrintperd, corrwtnpart) global(stddev, numperiod);

churnrate=numrepperd/clustsz;
corrclustperd=(clustsz*corrwtnperd*corrintperd+(1-churnrate)*(1-corrwtnperd)*corrwtnpart)/(1+(clustsz-1)*corrwtnperd);

f=999;
f=4*stddev**2/(numperiod*numclust*clustsz)*(1+(clustsz-1)*corrwtnperd)*(1+(numperiod-1)*corrclustperd);

if f<0 | f=999 then print churnrate corrclustperd f;

return (f);
finish var_ete_pacrt;

start var_ete_dec_pacrt(decclustsz) global(corrpar, stddev, numperiod, numrepperd, costperclust, costperpart, costperoutc, budget);

corrwtnperd = corrpar[1];
corrintperd = corrpar[2];
corrwtnpart = corrpar[3];

print corrpar declustsz;
f=999;

if decclustsz>0 then
do;
churnrate = numrepperd/decclustsz;
corrclustperd=(decclustsz*corrwtnperd*corrintperd+(1-corrwtnperd)*corrwtnpart*(1-churnrate))/(1+(decclustsz-1)*corrwtnperd);
decnumclust=budget/(costperclust+costperpart*(1+(numperiod-1)*churnrate)*decclustsz+costperoutc*numperiod*decclustsz);
f=4*stddev**2/(numperiod*decnumclust*decclustsz)*(1+(decclustsz-1)*corrwtnperd)*(1+(numperiod-1)*corrclustperd);
end;

if f<0 | f=999 then print corrwtnperd corrintperd corrwtnpart decclustsz f;
*print corrwtnperd corrintperd corrwtnpart corrclustperd decclustsz churnrate decnumclust f;
return (f);
finish var_ete_dec_pacrt;

start funcre_pacrt(x) global(corrpar, numclust, clustsz, minnumclust, maxnumclust, minclustsz, maxclustsz, numperiod, numrepperd, costperclust, costperpart, costperoutc, budget);
corrpar=j(3, 1);
corrpar[1:3] = x; * copy x -> corrpar;

*print numclust clustsz;

corrwtnperd = x[1];
corrintperd = x[2];
corrwtnpart = x[3];

churncost=costperpart*(numperiod-1)*numrepperd;

con=j(3, 3, .);

bound1=(budget/maxnumclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);
bound2=(budget/numclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);

con[1, 1]=max(minclustsz, numrepperd, int(bound1));
con[2, 1]=min(maxclustsz, int(bound2));

con[3, 1]=1;
con[3, 2]=1;
if corrwtnperd^=0 then con[3, 3]=1-1/corrwtnperd;

print con;

tc=j(13, 1, .);
tc[6]=1e-8;

x0=10;
optn=j(1, 10, .);

optn[1]=0;

call nlpqn(rc, decclustsz, "var_ete_dec_pacrt", x0, optn) blc=con;

decnumclust=budget/(costperclust+(costperpart+costperoutc*numperiod)*decclustsz+churncost);
minvar=var_ete_dec_pacrt(decclustsz);
temp_var=var_ete_pacrt(decnumclust, decclustsz, numrepperd, corrwtnperd, corrintperd, corrwtnpart);
actdeccost=decnumclust*(costperclust+(costperpart+costperoutc*numperiod)*decclustsz+churncost);

re=999;
actcost=numclust*(costperclust+(costperpart+costperoutc*numperiod)*clustsz+churncost);
var=var_ete_pacrt(numclust, clustsz, numrepperd, corrwtnperd, corrintperd, corrwtnpart);
if var>0 & minvar>=0 then re = minvar/var;

if re=999 | re>1 then
do;
print numrepperd decclustsz decnumclust actdeccost minvar;
print corrwtnperd corrintperd corrwtnpart numclust clustsz var actcost re;
end;

return (re);
finish funcre_pacrt;

start mmd_pacrt_int(corrwtnperdmin, corrwtnperdmax, corrintperdmin, corrintperdmax, corrwtnpartmin, corrwtnpartmax)
global(numclust, clustsz, costperclust, costperpart, costperoutc, budget, churncost,
numseq_lcrt, numperiod, numrepperd, clustszstart, minnumclust, maxnumclust, minclustsz, maxclustsz);

if mod(minnumclust, numseq_lcrt)=0 then numcluststart=minnumclust;
else numcluststart=(int(minnumclust/numseq_lcrt)+1)*numseq_lcrt;

bound1=(budget/maxnumclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);
clustszstart=max(minclustsz, numrepperd, int(bound1));

reminmax=0;

do numclust=numcluststart to maxnumclust by numseq_lcrt;
bound2=(budget/numclust-costperclust-churncost)/(costperpart+costperoutc*numperiod);
clustszend=min(maxclustsz, int(bound2));

do clustsz=clustszstart to clustszend;
con=j(2, 3, .);

con[1, 1]=corrwtnperdmin;
con[2, 1]=corrwtnperdmax;
con[1, 2]=corrintperdmin;
con[2, 2]=corrintperdmax;
con[1, 3]=corrwtnpartmin;
con[2, 3]=corrwtnpartmax;

par=j(10, 1, .);
par[9, 1]=4;

x0=j(1, 3, .);
x0[1, 1]=(corrwtnperdmin+corrwtnperdmax)/2;
x0[1, 2]=(corrintperdmin+corrintperdmax)/2;
x0[1, 3]=(corrwtnpartmin+corrwtnpartmax)/2;

optn = {0};
call nlpqn(rc, xres, "funcre_pacrt", x0, optn) blc=con par=par;
remin=funcre_pacrt(xres);
*print numclust clustsz xres remin reminmax;

if remin>reminmax then
do;
mmdcorrwtnperd=xres[1];
mmdcorrintperd=xres[2];
mmdcorrwtnpart=xres[3];
mmdnumclust=numclust;
mmdclustsz=clustsz;
mmdre=remin;
mmdactcost=numclust*(costperclust+(costperpart+costperoutc*numperiod)*clustsz+churncost);
reminmax=remin;
*print reflag numclust clustsz actcost remin revar retau reomega reminmax;
end; *end remin;
end; *end clustsz;
end; *end numclust;

mmdchurnrate=numrepperd/mmdclustsz;
mmd=corrwtnperdmin||corrwtnperdmax||corrintperdmin||corrintperdmax||corrwtnpartmin||corrwtnpartmax||mmdcorrwtnperd||mmdcorrintperd||mmdcorrwtnpart||numseq_lcrt||numperiod||numrepperd||mmdnumclust||mmdclustsz||mmdchurnrate||mmdre||mmdactcost;

return (mmd);
finish mmd_pacrt_int;

I got the following message:

ERROR: Matrix declustsz has not been set to a value.

statement : PRINT at line 3003 column 1
traceback : module VAR_ETE_DEC_PACRT at line 3003 column 1
module FUNCRE_PACRT at line 3003 column 1
module MMD_PACRT_INT at line 3003 column 1

ERROR: Execution error as noted previously. (rc=1000)

operation : NLPQN at line 3003 column 1
operands : *LIT1105, x0, optn, con

*LIT1105 1 row 1 col (character, size 17)

var_ete_dec_pacrt

x0 1 row 1 col (numeric)

10

optn 1 row 10 cols (numeric)

0 . . . . . .
. . .

con 3 rows 3 cols (numeric)

2 . .
100 . .
1 1 -7

statement : CALL at line 3003 column 1
traceback : module FUNCRE_PACRT at line 3003 column 1
module MMD_PACRT_INT at line 3003 column 1

statement : CALL at line 3003 column 1
traceback : module FUNCRE_PACRT at line 3003 column 1
module MMD_PACRT_INT at line 3003 column 1

ERROR: Execution error as noted previously. (rc=1000)

operation : NLPQN at line 3003 column 1
operands : *LIT1148, x0, optn, con, , par

*LIT1148 1 row 1 col (character, size 12)

funcre_pacrt

x0 1 row 3 cols (numeric)

0.125 0.035 0.35

optn 1 row 1 col (numeric)

0

con 2 rows 3 cols (numeric)

0.05 0.02 0.2
0.2 0.05 0.5

par 10 rows 1 col (numeric)

.
.
.
.
.
.
.
.
4
.

statement : CALL at line 3003 column 1
traceback : module MMD_PACRT_INT at line 3003 column 1

statement : CALL at line 3003 column 1
traceback : module MMD_PACRT_INT at line 3003 column 1

[Rick_SAS]

Did you generate the error messages from the program that you attached? The line numbers do not look correct.  Maybe the program that you pasted is part of a macro and the error log is from calling the macro?

While debugging your code, I strongly recommend that you avoid embedding code in a macro.  Only after everything in your program is working should you think about embedding the code in a macro call.

The SAS Log tells you where the error is occurring, and which function was called when the error occurred. Please read about how to interpret SAS IML error messages. Here are some tips for debugging IML programs:
1. There are two types of errors. Parse-time errors include improper syntax and are easier to find and fix. Run-time errors take more effort because they often depend on the data being analyzed. 

2. The Log tells you that the error occurs in the PRINT statement at line 3003.  The Log also tells you WHY the error occurs: because you are trying to print the matrix 'declustsz', but it has not been set to a value.

3. The Traceback gives the calling sequence for the function calls.

• The error occurs in the function VAR_ETE_DEC_PACRT, which is defined at line 3003.
• Who called the VAR_ETE_DEC_PACRT function? It was the function FUNCRE_PACRT, which is defined at line ????. Your log says line 3003, but that can't be correct because the previous message was from line 3003.
  
• Who called the FUNCRE_PACRT function? It was the MMD_PACRT_INT, which is defined at line ????.

4. We know that the error is in the VAR_ETE_DEC_PACR function. The next section of the Log tells you which line fails. It is the NLPQN call on line ????.

The problem in the NLPQN call occurs because the objective function is failing. The best way to figure out why is to thoroughly test and debug the objective function BEFORE you try to call it from an optimizer. 

Good luck with your debugging. 

[Xing_]

Thank you for the useful information. Now I've revised my code, please see the attached file, still with 4 functions. My goal is to find the minimum of RE (in "funcre_pacrt") across the values of (corrwtnperd, corrintperd, corrwtnpart), where the corresponding ranges are (corrwtnperdmin, corrwtnperdmax),  (corrintperdmin, corrintperdmax), (corrwtnpartmin, corrwtnpartmax), under a pair value of numclust and clustsz.

As you see, for numclust=52 and clustsz=3, nlpqn shows the min of RE is 1 at corrwtnperd=0.125, corrintperd=0.11, corrwtnpart=0.35.  However, I found RE=0.9234 at corrwtnperd=0.05, corrintperd=0.02, corrwtnpart=0.3. nlpqn does not provide the correct RE value.  

[Rick_SAS · Accepted solution]

You can discover more information about your problem and the NLP optimization by turning on diagnostic printing. For example, turn on minimal printing by setting

        optn = {0, 2};      

When you run the program with this option, you see that the algorithm is taking ZERO steps from the initial guess. It reports that the gradient of the objective function is the zero vector at that location. 

You can confirm that fact by calling the NLPFDD function to compute numerical derivatives of the objective function at the initial guess:

x0 = {0.125000 0.110000 0.350000};
call nlpfdd(f, grad, Hess, "funcre_pacrt", x0);
print f, grad, Hess;

You obtain that the gradient and the Hessian are both exactly 0. Clearly, the objective function is not correct. You must debug the objective function before you attempt to optimize it.