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.     ... View more

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 ... View more

I did check the manual and tried several ways. None of them work.    tc=j(13, 1, .); *tc[6]=0.000000000001; tc[7]=10; ... View more

cluster size can be decimal in the objective function.    Is there a way to set up a tight termination criterion or any other solution?  ... View more

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 ... View more

Got it! Thank you so much!  ... View more

One more question about  optn= {1}; Does it mean to solve a maximization problem? It aims to minimize the function, the optn[1] value should be 0 instead of 1? When I replace 1 by 0, the error message is shown as below.    ERROR: NLPHQN call: The number of observations must be set. Use the first element of the OPT vector. ERROR: (execution) Matrices do not conform to the operation. ... View more

Thank you for the code. It really helps!   The log file shows either "NOTE: ABSGCONV convergence criterion satisfied." or "NOTE: All parameters are actively constrained. Optimization cannot proceed." Is something wrong with the second note? It sounds that no optimization procedure was done. Any input here? ... View more

Thank you for the response. Yes, h(x, y, z) is a 1-dimention function for each x and y. The example I presented is simple since both functions are monotonic functions. The real ones are very complicated actually ☹️  ... View more