Number of iterations in NLPNM suproutine

[RaniaMamdouh]

Hello all;

How can i increase the number of iterations in proc IML, NLPNMS subroutine?

thanks

[Rick_SAS]

The doc for the termination criteria is here: SAS/IML(R) 12.3 User's Guide

tc[1] specifies the maximum number of iterations in the optimization process. To get, say, 5000 iterations, use the following:

tc = j(10,1,.); /* set all to default (missing) */

tc[1] = 5000;  /* override default maximum iterations */

[RaniaMamdouh]

Thanks for your reply

but when i apply your solution also i find the same following error

ERROR: NMSIMP Optimization cannot be completed.

ERROR: NMSIMP needs more than 1000 iterations or 3000 function calls.

what can i do in this proplem?

[Hutch_sas]

Looks like you need to increase the second control parameter, the maximum number of function calls, i.e.

tc[2] = 10000;

[Rick_SAS]

That is a difficult question to answer, as many books and papers have been written on the topic. A few things to try:

1) Look at the return code to see why it is failing.

2) Look at the iteration history to see whether the optimization is jumping around, diverging to infinity, stuck in a valley/ridge, or some other situation.

3) Try a better initial guess.

4) Relax some of the convergence criteria. If the objective function is flat (see #1 and #2), you might not be able to converge by using the default values.

5) Try a different algorithm.

[RaniaMamdouh]

I tried these things, but i think that the problem can be in writing the routine. I will briefly describe my problem, i minimize the "cost" and i have 4 linear constrains on "x's" and three non-linear constrains on  (arl), and this is the last part of my program which i think that something wrong in it:

tc = j(10,1,.); /* set all to default (missing) */

tc[1] = 5000;  /* override default maximum iterations */

con = { 0.   0.   0.  1.,

        1.   .    .   20.};

x={0.1  0.7  3.5  2};

optn= j(1,11,.);optn[1]=0;optn[2]= 5;optn[10]=3;optn[11]=0;

CALL nlpnms(rc,xres,"Cost",x,optn,con,tc) nlc="arl";

quit;

please tell me if there is wrong.

thanks

[Rick_SAS]

I don't see any linear constraints for x.  All I see are ranges.  You need four additional rows for the CON matrix.  See SAS/IML(R) 12.3 User's Guide

[RaniaMamdouh]

yes, i mean that ranges for x's not linear constrains. Is there any other wrong?

[Hutch_sas]

You might try overriding the second termination criteria, set tc[2] = 10000 and see if that affects your results. it appears that  you are running up against a limit on objective function calls.

[RaniaMamdouh]

i'm already tried to set tc[2]=10000 but nothing improved. Is the place of "tc" in the CALL line right or not?

[Rick_SAS]

Can you show the ARL module for nonlinear contstraints?

[RaniaMamdouh]

yes sure; it's as follows:

start arl(x);

      bb0=j(3,1,0.);

      m=100;

      shift=0;

* L is ucl;

*x[1] is lamda, x[2] is L, x[3] is k;

* t is number of regions;

t=2*m+1;

* delta is width of each region;

delta=2*x[2]/t;

*now we are constructing the vector z;

*vector z represents different locations for the chart statistics;

*these are the midpoints of the t intervals;

z=-x[2]+.5*delta;

do i=2 to t;

k=-x[2]+(i-.5)*delta;

z=z//k;

end;

* now we are calculating the incontrol transitional probability matrix r;

        do i=1 to t;

        zz=0;

        do j=1 to t;

        tau1=z[j, ]+(delta/2)-z[i, ];

        tau2=z[j, ]-(delta/2)-z[i, ];

if tau1<-(x[1]*x[3]) then fun1=tau1-(1-x[1])*x[3];

if abs(tau1)<=(x[1]*x[3])then fun1=tau1/x[1];

if tau1>(x[1]*x[3]) then fun1=tau1+(1-x[1])*x[3];

if tau2<-(x[1]*x[3]) then fun2=tau2-(1-x[1])*x[3];

if abs(tau2)<=(x[1]*x[3])then fun2=tau2/x[1];

if tau2>(x[1]*x[3]) then fun2=tau2+(1-x[1])*x[3];

              kone=z[i,]+fun1;

              ktwo=z[i,]+fun2;

              pij=probnorm(kone-shift)-probnorm(ktwo-shift);

              zz=zz//pij;

              end;

       if i=1 then r=zz;

       else

       r=r||zz;

       end;

       v=j(t,1,0);

       v[((t+1)/2), ]=1;

       v=t(v);

r=r[2:t+1,];

* now matrix r is t*t with the row of all zeros removed;

r=t(r);

*j(t,1,1) is tx1 vector of the value 1;

arl=inv((I(t)-r))*j(t,1,1);

arl0=v*arl;

bb0[1]=arl0-500;

m=100;

shift=1;

* L is ucl;

*x[1] is lamda, x[2] is L, x[3] is k;

* t is number of regions;

t=2*m+1;

* delta is width of each region;

delta=2*x[2]/t;

*now we are constructing the vector z;

*vector z represents different locations for the chart statistics;

*these are the midpoints of the t intervals;

z=-x[2]+.5*delta;

do i=2 to t;

k=-x[2]+(i-.5)*delta;

z=z//k;

end;

* now we are calculating the incontrol transitional probability matrix r;

        do i=1 to t;

        zz=0;

        do j=1 to t;

        tau1=z[j, ]+(delta/2)-z[i, ];

        tau2=z[j, ]-(delta/2)-z[i, ];

if tau1<-(x[1]*x[3]) then fun1=tau1-(1-x[1])*x[3];

if abs(tau1)<=(x[1]*x[3])then fun1=tau1/x[1];

if tau1>(x[1]*x[3]) then fun1=tau1+(1-x[1])*x[3];

if tau2<-(x[1]*x[3]) then fun2=tau2-(1-x[1])*x[3];

if abs(tau2)<=(x[1]*x[3])then fun2=tau2/x[1];

if tau2>(x[1]*x[3]) then fun2=tau2+(1-x[1])*x[3];

              kone=z[i,]+fun1;

              ktwo=z[i,]+fun2;

              pij=probnorm(kone-shift)-probnorm(ktwo-shift);

              zz=zz//pij;

              end;

       if i=1 then r=zz;

       else

       r=r||zz;

       end;

       v=j(t,1,0);

       v[((t+1)/2), ]=1;

       v=t(v);

r=r[2:t+1,];

* now matrix r is t*t with the row of all zeros removed;

r=t(r);

*j(t,1,1) is tx1 vector of the value 1;

arl=inv((I(t)-r))*j(t,1,1);

arl1=v*arl;

      bb0[2]=12-arl1;

m=100;

shift=0;

* L is ucl;

*x[1] is lamda, x[2] is L, x[3] is k, x[4] is n;

* t is number of regions;

t=2*m+1;

* delta is width of each region;

delta=2*x[2]/t;

*now we are constructing the vector z;

*vector z represents different locations for the chart statistics;

*these are the midpoints of the t intervals;

z=-x[2]+.5*delta;

do i=2 to t;

k=-x[2]+(i-.5)*delta;

z=z//k;

end;

* now we are calculating the incontrol transitional probability matrix r;

        do i=1 to t;

        zz=0;

        do j=1 to t;

        tau1=z[j, ]+(delta/2)-z[i, ];

        tau2=z[j, ]-(delta/2)-z[i, ];

if tau1<-(x[1]*x[3]) then fun1=tau1-(1-x[1])*x[3];

if abs(tau1)<=(x[1]*x[3])then fun1=tau1/x[1];

if tau1>(x[1]*x[3]) then fun1=tau1+(1-x[1])*x[3];

if tau2<-(x[1]*x[3]) then fun2=tau2-(1-x[1])*x[3];

if abs(tau2)<=(x[1]*x[3])then fun2=tau2/x[1];

if tau2>(x[1]*x[3]) then fun2=tau2+(1-x[1])*x[3];

              kone=z[i,]+fun1;

              ktwo=z[i,]+fun2;

              pij=probnorm(kone-shift)-probnorm(ktwo-shift);

              zz=zz//pij;

              end;

       if i=1 then r=zz;

       else

       r=r||zz;

       end;

       v=j(t,1,0);

       v[((t+1)/2), ]=1;

       v=t(v);

r=r[2:t+1,];

* now matrix r is t*t with the row of all zeros removed;

r=t(r);

*j(t,1,1) is tx1 vector of the value 1;

arl=inv((I(t)-r))*j(t,1,1);

arl0=v*arl;

m=100;

shift=4;

* L is ucl;

*x[1] is lamda, x[2] is L, x[3] is k;

* t is number of regions;

t=2*m+1;

* delta is width of each region;

delta=2*x[2]/t;

*now we are constructing the vector z;

*vector z represents different locations for the chart statistics;

*these are the midpoints of the t intervals;

z=-x[2]+.5*delta;

do i=2 to t;

k=-x[2]+(i-.5)*delta;

z=z//k;

end;

* now we are calculating the incontrol transitional probability matrix r;

        do i=1 to t;

        zz=0;

        do j=1 to t;

        tau1=z[j, ]+(delta/2)-z[i, ];

        tau2=z[j, ]-(delta/2)-z[i, ];

if tau1<-(x[1]*x[3]) then fun1=tau1-(1-x[1])*x[3];

if abs(tau1)<=(x[1]*x[3])then fun1=tau1/x[1];

if tau1>(x[1]*x[3]) then fun1=tau1+(1-x[1])*x[3];

if tau2<-(x[1]*x[3]) then fun2=tau2-(1-x[1])*x[3];

if abs(tau2)<=(x[1]*x[3])then fun2=tau2/x[1];

if tau2>(x[1]*x[3]) then fun2=tau2+(1-x[1])*x[3];

              kone=z[i,]+fun1;

              ktwo=z[i,]+fun2;

              pij=probnorm(kone-shift)-probnorm(ktwo-shift);

              zz=zz//pij;

              end;

       if i=1 then r=zz;

       else

       r=r||zz;

       end;

       v=j(t,1,0);

       v[((t+1)/2), ]=1;

       v=t(v);

r=r[2:t+1,];

* now matrix r is t*t with the row of all zeros removed;

r=t(r);

*j(t,1,1) is tx1 vector of the value 1;

arl=inv((I(t)-r))*j(t,1,1);

arl1=v*arl;

c0=10;

c1=100;

gam1=1;

gam2=1;

O=0.01;

W=25;

a=0.5;

b=0.1;

F=50;

E=0.05;

T0=0;

T1=2;

T2=0;

h4=(x[4]*E)+(gam1*T1)+(gam2*(T2**2));

h1=((arl1-0.5)*(O*(F+(c1*T0*(gam1-1)))-(2*arl0*(c0+O*(((arl1-0.5)*a)+((arl1-0.5)*b*x[4])+W)+(c1*(-1+h4*O-O*((x[4]*E)+T1+T2)))))))/(2*O*arl0);

h2=((-2*(arl1-0.5)*(F+(c1*T0*(gam1-1))+(arl0*(a+(b*x[4]))*(1+(h4*O))))))/(O*arl0);

h3=(-1/(2*(O**2)*arl0))*((2*F)+(2*c0*T0*(gam1-1))-(a*T0*O)-(2*(arl1-0.5)*a*T0*O)-(2*c1*h4*T0*O)-(b*x[4]*T0*O)-(2*(arl1-0.5)*b*x[4]*T0*O)-(2*T0*W*O)+(2*F*O*((x[4]*E)+T1+T2))+(a*T0*gam1*O)+(2*(arl1-0.5)*a*T0*gam1*O)+(2*c1*h4*T0*gam1*O)+(b*x[4]*T0*gam1*O)+(2*(arl1-0.5)*b*x[4]*T0*gam1*O)+(2*T0*W*gam1*O)-(a*h4*T0*(O**2))-(b*x[4]*h4*T0*(O**2))+(a*h4*T0*gam1*(O**2))+(a*x[4]*h4*T0*gam1*(O**2))+(2*arl0*(a+b*x[4])*(1+h4*O)*(1+O*((x[4]*E)+T1+T2))));

h=(-h2+SQRT((h2**2)-(4*h1*h3)))/(2*h1);

s=(1/(O*h))-0.5;

r1=(c0/O)+(c1*((x[4]*E)+(h*(arl1-0.5))+(gam1*T1)+(gam2*T2))+((s*F)/arl0)+W);

r2=((a+(b*x[4]))/h)*((1/O)+(x[4]*E)+(h*(arl1-0.5))+(gam1*T1)+(gam2*T2));

r3=(1/O)+((1-gam1)*s*T0)/arl0+(x[4]*E)+(h*(arl1-0.5))+T1+T2;

C2=(r1+r2)/r3;

print C2;

bb0[3]=(1.005*13.355)-C2;

    return (bb0);

finish arl;

[Rick_SAS]

Ugh. Too complicated. I don't see how to verify that the ARL module is correct.

Can you solve the unconstrained problem? What about with only the first or only the first 2 constraints?

[RaniaMamdouh]

if i put only the third constrain i can face the same problem, i.e. the problem isn't in the number of constrains?

[Rick_SAS]

Is it possible that the problem is in the (very complicated) nonlinear constraints? Can you show that there is a feasible solution that satisfies the constraints? If not, you will encounter this kind of non-convergence.

To give a simple example, the problem min F(x,y) s.t. y> -x^2 and y>x^2 +1 has no solution because no (x,y) pair is feasible.