BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
dpullarkat
Calcite | Level 5

Hi,

 

I have already looked at the other threads on the SAS community dealing with similar problems; and I have already tried those solutions

  • one which suggested for the objective to be multiplied with 1.00e-6, so that the scale of the objective would be smaller
  • another dealing with RELOPTOL and FDCENTRAL

Both have not worked.

 

I am sure that the logic is correct, since something very similar is being currently run in excel using the Solver add-in.

 

The code is attached below (alongwith the datasets). Would be very grateful if someone could help!

 

Thanks!

 

/*read in units, dot and initial compliance data*/

data drug_consumption;
length drug $ 14;
input drug $ dailydose compliance;
datalines;
drug1	 0.14 	 0.71 
drug2	 5.11 	 0.61 
drug3	 6.00 	 0.50 
drug4	 12.00 	 0.50 
;
run;

data regimen_dot;
length regimen $ 16;
input regimen $ dot;
datalines;
reg1	37
reg2	42
reg3	41
;
run;

data regimen_drug;
length regimen $16;
length drug $14;
input regimen drug weight;
datalines;
reg1	drug1	1
reg2	drug1	1
reg3	drug1	1
reg1	drug2	1
reg2	drug2	1
reg3	drug2	1
reg3	drug3	0.285714286
reg2	drug4	0.292682927
reg3	drug4	0
reg2	drug3	0
reg1	drug4	0
reg1	drug3	0
;
run;

data units_data;
length drug $14.;
input drug $ month $ units;
datalines;
drug1	12-Jan	15729
drug1	12-Feb	15019
drug1	12-Mar	17511
drug1	12-Apr	14647
drug1	12-May	16263
drug1	12-Jun	13807
drug1	12-Jul	15235
drug1	12-Aug	14350
drug1	12-Sep	13301
drug1	12-Oct	15822
drug1	12-Nov	16321
drug1	12-Dec	14393
drug2	12-Jan	506782
drug2	12-Feb	476800
drug2	12-Mar	555050
drug2	12-Apr	452674
drug2	12-May	526190
drug2	12-Jun	443008
drug2	12-Jul	488768
drug2	12-Aug	460306
drug2	12-Sep	420596
drug2	12-Oct	510478
drug2	12-Nov	515308
drug2	12-Dec	458582
drug3	12-Jan	8064
drug3	12-Feb	7308
drug3	12-Mar	9240
drug3	12-Apr	12726
drug3	12-May	19992
drug3	12-Jun	28728
drug3	12-Jul	30954
drug3	12-Aug	37800
drug3	12-Sep	44016
drug3	12-Oct	88032
drug3	12-Nov	87528
drug3	12-Dec	54558
drug4	12-Jan	17808
drug4	12-Feb	15456
drug4	12-Mar	20496
drug4	12-Apr	16800
drug4	12-May	24864
drug4	12-Jun	18144
drug4	12-Jul	27552
drug4	12-Aug	38640
drug4	12-Sep	42000
drug4	12-Oct	58464
drug4	12-Nov	79296
drug4	12-Dec	90048
;
run;

data month;
input month $ ;
datalines;
12-Jan
12-Feb
12-Mar
12-Apr
12-May
12-Jun
12-Jul
12-Aug
12-Sep
12-Oct
12-Nov
12-Dec
;
run;


PROC OPTMODEL;

   SET <STRING> DRUG;
   SET <STRING> REGIMEN;
   SET <STRING> MONTH;

   read data month into month=[month];

   NUMBER dailydose{DRUG};
   READ DATA DRUG_consumption INTO DRUG=[DRUG] dailydose;

   NUMBER DOT{regimen};
   READ DATA REGIMEN_dot INTO regimen=[regimen] dot;    
   PRINT DOT;

   NUMBER UNITS{DRUG, MONTH};
   READ DATA units_data INTO [DRUG MONTH] UNITS;
   PRINT UNITS;

   number weight{regimen, drug};
   read data regimen_drug into [regimen drug] weight;
   print weight;
  
   /*DEFINE THE PARAMETERS*/

   var ps {REGIMEN, MONTH}  init 0.35 >=0.01 <=1,
   annualpatients init 10000 <=9400 >=4558, 
   compliance {drug} init 0.1 >=0.3 <=1;
   
   
   num su_add{ d in drug, m in month}=units[d,m]/dailydose[d];

   impvar monthly_patient_share{m in month}=sum{r in regimen}ps[r,m];
   impvar regimen_tot_ps{r in regimen}=sum{m in month, r in regimen} ps[r,m];
   impvar annual_ps{r in regimen}=regimen_tot_ps[r]/12;
   impvar avg_dot_1=sum(sum{m in month, r in regimen} (ps[r,m]*dot[r]));
   impvar avg_dot=avg_dot_1/12;
   impvar avg_dot_months=avg_dot/(30/7);

/*   annual total treated patients*/
   impvar pat_tot=annualpatients*avg_dot_months;
   impvar total_patients=sum{m in month, r in regimen} (ps[r,m]*annualpatients/(avg_dot/(30/7)));
   num total_units=sum{m in month,d in drug} su_add[d,m];
   num monthly_units{m in month}=sum{d in drug} su_add[d,m];
   num ratio{m in month}=monthly_units[m]/total_units;
   impvar monthlypatients{m in month}=pat_tot*ratio[m];
   
/*   calculate total drug->regimen shares*/
   impvar drug_share{d in drug, m in month} = sum{m in month, r in regimen}ps[r,m]*weight[r,d];
/*   calculated units*/
    impvar calc_units{d in drug, m in month}=monthlypatients[m]*dailydose[d]*compliance[d]*30*drug_share[d,m];
/*	deviation*/
    impvar deviation{d in drug, m in month}=abs(units[D,M]-calc_units[d,m])/calc_units[d,m];

   /*SET THE CONSTRAINTS*/
   constraint  MONTHLY_ps {M IN MONTH}:  monthly_patient_share[m]=1;
   constraint  patient_match: total_patients=annualpatients;   
   constraint  calc_units_dev{d in drug, m in month}: deviation[d,m]<=0.01;


   /*SET THE OBJECTIVE*/
   min error = sum{M IN MONTH, d in drug}( ((units[D,M]-calc_units[d,m])/(dailydose[d]))**2 )/10000000000000000;
  
 expand;
 solve with nlpc/  tech=quanew;
 print deviation;
 print ps;
 print compliance;
 print monthly_patient_share;
 print annualpatients;
 print total_patients;
 print calc_units;

 quit;
1 ACCEPTED SOLUTION

Accepted Solutions
RobPratt
SAS Super FREQ

I have a few more questions and suggestions.

 

Why do you initialize some of the variables to values outside their bounds?

 

The NLPC solver was retired (no longer documented or maintained) in SAS 9.3.  You should instead use the NLP solver:

http://support.sas.com/documentation/cdl/en/ormpug/68156/HTML/default/viewer.htm#ormpug_nlpsolver_to...

 

It is better to avoid using the absolute value function, which is non-differentiable.  Instead, you can rewrite those statements as:

   impvar deviation{d in drug, m in month}=(units[D,M]-calc_units[d,m])/calc_units[d,m];
   constraint calc_units_dev{d in drug, m in month}: -0.1 <= deviation[d,m] <= 0.1;

 

For sum-of-squares objectives, it is sometimes better to introduce explicit variables so that you can express the objective as a sum of squares of the new variables, yielding a sparse (diagonal) Hessian:

   var error_m_d {M IN MONTH, d in drug};
   con error_con {M IN MONTH, d in drug}: error_m_d[m,d] = (units[D,M]-calc_units[d,m])/dailydose[d];
   min error = sum{M IN MONTH, d in drug} error_m_d[m,d]**2;

View solution in original post

5 REPLIES 5
RobPratt
SAS Super FREQ

First, let's make sure the model is correct.  I see a couple of issues.

 

The "r in regimen" should not appear on both sides of this declaration:

impvar regimen_tot_ps{r in regimen}=sum{m in month, r in regimen} ps[r,m];

 

The "m in month" should not appear on both sides of this declaration:

 

impvar drug_share{d in drug, m in month} = sum{m in month, r in regimen}ps[r,m]*weight[r,d];

 

Is there a reason for the outer SUM here?

impvar avg_dot_1=sum(sum{m in month, r in regimen} (ps[r,m]*dot[r]));

After correcting these, do the results from EXPAND look like what you expect?

 

As another check on the correctness of the model, you should probably FIX the variables to their values from the Excel solution and then call the solver.

 

dpullarkat
Calcite | Level 5

Hi Robb,

 

Thanks for looking into the issue. I've implemented the changes you've suggested and also tried fixing the variable to those from the excel solver. All values of the imputed variables and the objective match in that case, so the logic of the optmodel is fine.

 

However, what is confusing me is why the solver is not trying the run more iterations. When I check the output (without fixing the values), the # of iterations is 0.

 

I've attached the updated code.

 

Thanks,

Deepthi

/*read in units, dot and initial compliance data*/

data drug_consumption;
length drug $ 14;
input drug $ dailydose compliance;
datalines;
drug1	 0.14 	 0.71 
drug2	 5.11 	 0.61 
drug3	 6.00 	 0.50 
drug4	 12.00 	 0.50 
;
run;

data regimen_dot;
length regimen $ 16;
input regimen $ dot;
datalines;
reg1	37
reg2	42
reg3	41
;
run;

data regimen_drug;
length regimen $16;
length drug $14;
input regimen drug weight;
datalines;
reg1	drug1	1
reg2	drug1	1
reg3	drug1	1
reg1	drug2	1
reg2	drug2	1
reg3	drug2	1
reg3	drug3	0.285714286
reg2	drug4	0.292682927
reg3	drug4	0
reg2	drug3	0
reg1	drug4	0
reg1	drug3	0
;
run;

data units_data;
length drug $14.;
input drug $ month $ units;
datalines;
drug1	12-Jan	15729
drug1	12-Feb	15019
drug1	12-Mar	17511
drug1	12-Apr	14647
drug1	12-May	16263
drug1	12-Jun	13807
drug1	12-Jul	15235
drug1	12-Aug	14350
drug1	12-Sep	13301
drug1	12-Oct	15822
drug1	12-Nov	16321
drug1	12-Dec	14393
drug2	12-Jan	506782
drug2	12-Feb	476800
drug2	12-Mar	555050
drug2	12-Apr	452674
drug2	12-May	526190
drug2	12-Jun	443008
drug2	12-Jul	488768
drug2	12-Aug	460306
drug2	12-Sep	420596
drug2	12-Oct	510478
drug2	12-Nov	515308
drug2	12-Dec	458582
drug3	12-Jan	8064
drug3	12-Feb	7308
drug3	12-Mar	9240
drug3	12-Apr	12726
drug3	12-May	19992
drug3	12-Jun	28728
drug3	12-Jul	30954
drug3	12-Aug	37800
drug3	12-Sep	44016
drug3	12-Oct	88032
drug3	12-Nov	87528
drug3	12-Dec	54558
drug4	12-Jan	17808
drug4	12-Feb	15456
drug4	12-Mar	20496
drug4	12-Apr	16800
drug4	12-May	24864
drug4	12-Jun	18144
drug4	12-Jul	27552
drug4	12-Aug	38640
drug4	12-Sep	42000
drug4	12-Oct	58464
drug4	12-Nov	79296
drug4	12-Dec	90048
;
run;

data month;
input month $ ;
datalines;
12-Jan
12-Feb
12-Mar
12-Apr
12-May
12-Jun
12-Jul
12-Aug
12-Sep
12-Oct
12-Nov
12-Dec
;
run;


PROC OPTMODEL;

   SET <STRING> DRUG;
   SET <STRING> REGIMEN;
   SET <STRING> MONTH;

   read data month into month=[month];

   NUMBER dailydose{DRUG};
   READ DATA DRUG_consumption INTO DRUG=[DRUG] dailydose;

   NUMBER DOT{regimen};
   READ DATA REGIMEN_dot INTO regimen=[regimen] dot;    
   PRINT DOT;

   NUMBER UNITS{DRUG, MONTH};
   READ DATA units_data INTO [DRUG MONTH] UNITS;
   PRINT UNITS;

   number weight{regimen, drug};
   read data regimen_drug into [regimen drug] weight;
   print weight;
  
   /*DEFINE THE PARAMETERS*/

   var ps {REGIMEN, MONTH}  init 0.35 >=0.01 <=1,
   annualpatients init 10000 <=9400 >=4558, 
   compliance {drug} init 0.1 >=0.3 <=1;
   
   
   num su_add{ d in drug, m in month}=units[d,m]/dailydose[d];

 impvar monthly_patient_share{m in month}=sum{r in regimen}ps[r,m];
   impvar regimen_tot_ps{r in regimen}=sum{m in month} ps[r,m];
   impvar annual_ps{r in regimen}=regimen_tot_ps[r]/12;
   impvar avg_dot_1 {r in regimen}=annual_ps[r]*dot[r];
   impvar avg_dot_2=sum{r in regimen} avg_dot_1[r];
   impvar avg_dot=avg_dot_2;
   impvar avg_dot_months=avg_dot_2/(30/7);

/*   annual total treated patients*/
   impvar pat_tot=annualpatients*avg_dot_months;
   impvar total_patients=sum{m in month, r in regimen} (ps[r,m]*annualpatients/(avg_dot/(30/7)));
   num total_units=sum{m in month,d in drug} su_add[d,m];
   num monthly_units{m in month}=sum{d in drug} su_add[d,m];
   num ratio{m in month}=monthly_units[m]/total_units;
   impvar monthlypatients{m in month}=pat_tot*ratio[m];
   
/*   calculate total drug->regimen shares*/
   impvar drug_share{d in drug, m in month} = sum{r in regimen}ps[r,m]*weight[r,d];
/*   calculated units*/
    impvar calc_units{d in drug, m in month}=monthlypatients[m]*dailydose[d]*compliance[d]*30*drug_share[d,m];
/*	deviation*/
    impvar deviation{d in drug, m in month}=abs(units[D,M]-calc_units[d,m])/calc_units[d,m];

   /*SET THE CONSTRAINTS*/
   constraint  MONTHLY_ps {M IN MONTH}:  monthly_patient_share[m]=1;
   constraint  patient_match: total_patients=annualpatients;   
   constraint  calc_units_dev{d in drug, m in month}: deviation[d,m]<=0.1;

   /*SET THE OBJECTIVE*/
   min error = sum{M IN MONTH, d in drug}( ((units[D,M]-calc_units[d,m])/(dailydose[d]))**2 );
    

/*min error = sum{d in drug, m in month} (abs(units[D,M]-calc_units[d,m])/calc_units[d,m]);*/

 expand;
 solve with nlpc/  tech=quanew;
 print ps;
 quit;
RobPratt
SAS Super FREQ

I have a few more questions and suggestions.

 

Why do you initialize some of the variables to values outside their bounds?

 

The NLPC solver was retired (no longer documented or maintained) in SAS 9.3.  You should instead use the NLP solver:

http://support.sas.com/documentation/cdl/en/ormpug/68156/HTML/default/viewer.htm#ormpug_nlpsolver_to...

 

It is better to avoid using the absolute value function, which is non-differentiable.  Instead, you can rewrite those statements as:

   impvar deviation{d in drug, m in month}=(units[D,M]-calc_units[d,m])/calc_units[d,m];
   constraint calc_units_dev{d in drug, m in month}: -0.1 <= deviation[d,m] <= 0.1;

 

For sum-of-squares objectives, it is sometimes better to introduce explicit variables so that you can express the objective as a sum of squares of the new variables, yielding a sparse (diagonal) Hessian:

   var error_m_d {M IN MONTH, d in drug};
   con error_con {M IN MONTH, d in drug}: error_m_d[m,d] = (units[D,M]-calc_units[d,m])/dailydose[d];
   min error = sum{M IN MONTH, d in drug} error_m_d[m,d]**2;
dpullarkat
Calcite | Level 5

Robb,

 

Thanks a ton!

 

I was using SAS 9.2, which I think did not have nlp. I downloaded SAS 9.3 and tried the nlp solver and it worked! 

wezhou
SAS Employee

I suspect that this problem may be infeasible. 

 

For nonlinear optimization, determining if a problem is globally infeasible is equivalent to finding a global solution of a bound least square problem.  For such a reformulation, we can apply the local solvers as well as multi-start.  If the objective is 0, then the problem has feasible points.  If the objective is nonzero, then the problem may be infeasible.  In general we cannot guarantee that a feasible region does not exist somewhere, as the corresponding problem is NP-hard.

 

For example, I have created a new problem from your original sas file. First, I have converted your inequality constraints calc_units_dev{d in drug, m in month}: deviation[d,m]<=0.01 to the equality constraints: calc_units_dev{d in drug, m in month}: deviation[d,m] + newslack [d,m]=0.01 by introducing new slack variables newslack[d,m] >= 0. Your objective function has then been replaced by error = (total_patients-annualpatients)^2 + sum{M IN MONTH}(monthly_patient_share[m]-1)^2+ sum{d in drug, m in month}(newslack[d,m] + deviation[d,m] - 0.01)^2, and your all the constraints has been dropped out. Essentially, I have created a new least square problem with the box constraints to detect the feasibility of your original problem. If the optimal error is 0, then the original problem is feasible. The log file that I have got, with V9.3 is:

         Objective                           Optimality

Iter    Value       Infeasibility        Error

0       7138983   0                       40.89185820

1       4975949   0                       0.00500000

2       4612514   0.00514115      0.00514115

3      1218025    0                       0.00500000

4      1149047    0                       0.00005000

5      1065987    0                       0.00006865

6      1065901    0.13030618      0.13030618

7      1065767    0                       0.00005000

8      1068925    0                       0.00005000

9      1068107   0                       0.0000005021859

NOTE: Optimal.

NOTE: Objective = 1068107.19.

Thus, no a feasible point is found. The sas file named ReformulatedOrigianl for this new created problem has been attached.

 

I then removed your constraints  patient_match and calc_units_dev from your problem. Instead, only constraint MONTHLY_ps is considered. I have created another sas file called ReformulatedFeasible and it has been attached as well. The solver has then found an optimal solution with objective value being 0. That is, a feasible point has been found.

 

To further investigate if a feasible point can be found from different starting points, I ran the solver with MultiStart option on. It basically tries a nlp solver from many different starting points. The command is:

solve with nlp/MS;

 

But, no a feasible point has been found. So I suspect that your original problem may be infeasible.

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 5 replies
  • 1233 views
  • 0 likes
  • 3 in conversation