Hi,
I have already looked at the other threads on the SAS community dealing with similar problems; and I have already tried those solutions
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;
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:
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;
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.
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;
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:
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;
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!
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.
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
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.