Proc Optmodel Problem Scaling

[svo]

Hi-

I'm trying to use proc optmodel and scale beyond 1 set of parameters for optimization (i.e. extend a 1 case example to a dataset). In the first instance I've written the following to find s such that the difference between the variables newrate and r are minimized:

/* Proc Optmodel */

proc optmodel;
number b=.33267, a1=.4361836, a2=-.1201676, a3=.937298, ey=198, yg=168.3, r=.038;
var s init ey, z, t, p, el;

min newrate = r;

con z=(1/sqrt(2*3.14))*exp(-.5*((ey-yg)/s)**2);
con t=1/(1+b*((ey-yg)/s));
con p=z*(a1*t+a2*t**2+a3*t**3);
con el=p*(yg-ey)+z*s;
con newrate=el/yg;

solve;
print s;

However I am having trouble extending this to the case where I have multiple observations for parameters ey, yg, r, and I am picking S for each observation.  Here's as far as I've gotten:

data parms;

input name $ ey yg b a1 a2 a3 r;

datalines;

COL1 100 80 .33267 .4361836 -.120167 .937298 .038

COL2 150 130 .33267 .4361836 -.120167 .937298 .038

COL3 100 75 .33267 .4361836 -.120167 .937298 .038

run;

proc optmodel;

set <string> parameters;

number ey{parameters};

number yg{parameters};

number b{parameters};

number a1{parameters};

number a2{parameters};

number a3{parameters};

number r{parameters};

read data parms into parameters=[name] ey=ey yg=yg b=b a1=a1 a2=a2 a3=a3 r=r;

var s{parameters} init 10;

var z{parameters};

var t{parameters};

var pr{parameters};

var el{parameters};

minimize newrate = r;

con z{p in parameters} = (1/sqrt(2*3.14))*exp(-.5*((ey-yg)/s)**2);

con t{p in parameters} = 1/(1+b*((ey-yg)/s));

con pr{p in parameters} =z*(a1*t+a2*t**2+a3*t**3);

con el{p in parameters}=p*(yg-ey)+z*s;

con newrate=el/yg;

solve;

print s;

Any help would be much appreciated!

Thanks,
SVO

[RobPratt]

To clarify, do you want to solve a separate optimization problem for each observation in the parms data set?

Also, your model does not minimize the difference between newrate and r.  Instead, it forces newrate to be equal to r.  Your objective function is a constant, and the following constraint would achieve the same purpose:

con r=el/yg;

[svo]

Hi Rob-

Yes, that's exactly what I'm trying to do.  Solve for the value of s such that new rate =/~ rate for each observation where ey and yg vary (r could also vary but not in my example code).  

Appreciate any thoughts on your end.

[RobPratt · Accepted solution]

Here is a way to solve one problem per observation and store the resulting values of s:

proc optmodel;
   set <string> parameters;
    
   number ey{parameters};
   number yg{parameters};
   number b{parameters};
   number a1{parameters};
   number a2{parameters};
   number a3{parameters};
   number r{parameters};
    
   read data parms into parameters=[name] ey=ey yg=yg b=b a1=a1 a2=a2 a3=a3 r=r;
    
   var s init 10;
   var z;
   var t;
   var pr;
   var el;

   str p; 
   minimize newrate = r[p];
    
   con z = (1/sqrt(2*3.14))*exp(-.5*((ey[p]-yg[p])/s)**2);
   con t = 1/(1+b[p]*((ey[p]-yg[p])/s));
   con pr =z*(a1[p]*t+a2[p]*t**2+a3[p]*t**3);
   con el=pr*(yg[p]-ey[p])+z*s;
   con newrate=el/yg[p];

   num s_sol {parameters};
   do p = parameters; 
      solve;
      print s;
      s_sol[p] = s;
   end;
   print s_sol;
quit;

You can instead solve these independent problems in parallel by replacing the DO loop with a COFOR loop:

   cofor {pp in parameters} do;
      p = pp;
      solve;
      print s;
      s_sol[p] = s;
   end;

A similar idea is illustrated in this documentation example.

[svo]

Thanks Rob- appreciate it.  This is exactly what I was trying to achieve.  

[svo]

Rob-

I spoke to fast.  I'm not sure why (maybe I am misunderstanding something) but I am seeing some weird behavior when I extend this to a large group of observations.  For instance:

/* Proc Optmodel Example One*/

data parms;
input name $ ey yg b a1 a2 a3 r;
datalines;
ID1 100 80 .33267 .4361836 -.120167 .937298 .038
ID2 150 130 .33267 .4361836 -.120167 .937298 .038
ID3 100 75 .33267 .4361836 -.120167 .937298 .038
ID4 198 168.63 .33267 .4361836 -.120167 .937298 .038
ID5 1053 684.45 .33267 .4361836 -.120167 .937298 .0301141374
run;

proc optmodel;
set <string> parameters;

number ey{parameters};
number yg{parameters};
number b{parameters};
number a1{parameters};
number a2{parameters};
number a3{parameters};
number r{parameters};

read data parms into parameters=[name] ey=ey yg=yg b=b a1=a1 a2=a2 a3=a3 r=r;


var s init 10;
var z;
var t;
var pr;
var el;

str p;
minimize newrate = r[p];

con z = (1/sqrt(2*3.14))*exp(-.5*((ey[p]-yg[p])/s)**2);
con t = 1/(1+b[p]*((ey[p]-yg[p])/s));
con pr =z*(a1[p]*t+a2[p]*t**2+a3[p]*t**3);
con el=pr*(yg[p]-ey[p])+z*s;
con newrate=el/yg[p];

/* Do Loop Version */

num s_sol_do {parameters};
do p = parameters;
solve;
s_sol_do[p] = s;
end;
create data solution_do from [name] s_sol_do;

/* Parallel Example */
num s_sol_cofor {parameters};
cofor {pp in parameters} do;
p = pp;
solve;
s_sol_cofor[p] = s;
create data solution_cofor from [name] s_sol_cofor;
end;
quit;

proc print data=solution_cofor; run;

Obs name cofor

1 ID1 25.118
2 ID2 31.258
3 ID3 28.001
4 ID4 43.266
5 ID5 323.855


proc print data=solution_do; run;

Obs name s_sol_do

1 ID1 25.118
2 ID2 31.258
3 ID3 28.001
4 ID4 43.266
5 ID5 123.713

So ID5 has a different answer in the do version vs cofor version.  Using a different optimization package in R and brute force method in SAS agree that the answer should be around 323.8 (the cofor version).  This is also confirmed in the one case verion below:

/* Proc Optmodel */

proc optmodel;
number b=.33267, a1=.4361836, a2=-.1201676, a3=.937298, ey=1053, yg=684.45, r=.0301141374;
var s init ey, z, t, p, el;

min newrate = r;

con z=(1/sqrt(2*3.14))*exp(-.5*((ey-yg)/s)**2);
con t=1/(1+b*((ey-yg)/s));
con p=z*(a1*t+a2*t**2+a3*t**3);
con el=p*(yg-ey)+z*s;
con newrate=el/yg;

solve;
print s;

/* OUTPUT */

Solver NLP
Algorithm Interior Point
Objective Function newrate
Solution Status Optimal
Objective Value 0.0301141374

Optimality Error 8.4013373E-8
Infeasibility 8.4013373E-8

Iterations 14
Presolve Time 0.00
Solution Time 0.04

s

323.85

However, suppose I add another two lines of input data to the parms dataset:

data parms;
input name $ ey yg b a1 a2 a3 r;
datalines;
ID1 100 80 .33267 .4361836 -.120167 .937298 .038
ID2 150 130 .33267 .4361836 -.120167 .937298 .038
ID3 100 75 .33267 .4361836 -.120167 .937298 .038
ID4 198 168.63 .33267 .4361836 -.120167 .937298 .038
ID5 1053 684.45 .33267 .4361836 -.120167 .937298 .0301141374
ID6 1180 826 .33267 .4361836 -.120167 .937298 .1345990522
ID7 100 80 .33267 .4361836 -.120167 .937298 .1345990522
run;

proc optmodel;
set <string> parameters;

number ey{parameters};
number yg{parameters};
number b{parameters};
number a1{parameters};
number a2{parameters};
number a3{parameters};
number r{parameters};

read data parms into parameters=[name] ey=ey yg=yg b=b a1=a1 a2=a2 a3=a3 r=r;


var s init 10;
var z;
var t;
var pr;
var el;

str p;
minimize newrate = r[p];

con z = (1/sqrt(2*3.14))*exp(-.5*((ey[p]-yg[p])/s)**2);
con t = 1/(1+b[p]*((ey[p]-yg[p])/s));
con pr =z*(a1[p]*t+a2[p]*t**2+a3[p]*t**3);
con el=pr*(yg[p]-ey[p])+z*s;
con newrate=el/yg[p];

/* Do Loop Version */

num s_sol_do {parameters};
do p = parameters;
solve;
s_sol_do[p] = s;
end;
create data solution_do from [name] s_sol_do;

/* Parallel Example */
num s_sol_cofor {parameters};
cofor {pp in parameters} do;
p = pp;
solve;
s_sol_cofor[p] = s;
create data solution_cofor from [name] s_sol_cofor;
end;
quit;

proc print data=solution_cofor; run;

Obs name cofor

1 ID1 -5.67537
2 ID2 -5.75878
3 ID3 -6.55722
4 ID4 -7.71885
5 ID5 -5.71067
6 ID6 -5.71067
7 ID7 -5.88630

proc print data=solution_do; run;

Obs name s_sol_do

1 ID1 25.118
2 ID2 31.258
3 ID3 28.001
4 ID4 43.266
5 ID5 123.713
6 ID6 624.578
7 ID7 -5.711

So some very weird answers and I'm not sure what I've done to send things haywire like this.  Would appreciate any thoughts.

Thanks,

[RobPratt]

I'm not sure anything is wrong here.  Your objective function is a constant, so all feasible solutions are equally desirable.  Different starting solutions can lead to different optimal solutions.  When you solve serially (with the DO loop), each solve starts from the previous solution.  When you run in parallel (with COFOR), the solves are independent.

If you have some way to quantify that one feasible solution is preferable to another, you can make that your objective instead.  For example, if you prefer larger values of s, you can make the following changes:

/*minimize newrate = r[p];*/
maximize myobj = s;

/*con newrate=el/yg[p];*/
con r[p]=el/yg[p];

[svo]

My misunderstanding, I needed to impose two additional constraints including:

con pr >= 0

con s >= 0

to achieve the desired result.  Thanks again.  

[RobPratt]

OK.  By the way, you can impose those bounds in the variable declarations, without having to introduce explicit constraints, as follows:

var pr >= 0;
var s >= 0;

[svo]

Got it.  Thanks