SAS Optimization, and SAS Simulation Studio

turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

Topic Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-23-2016 01:20 AM

I have a Proc Optmodel code which runs (on SAS 9.3), but the solved values of the variables which I have are beyond the upper and lower bounds I have specified in my var declaration.

Any idea why this could be happening and what I can do to fix it?

As you can see from the code below, ps can been bounded between 0 and 1. But in the solved output, the values frequently go negative or above 1. I am having a similar issue with compliance, which also goes above 1 or negative.

```
PROC OPTMODEL;
SET <STRING> DRUG;
SET <STRING> REGIMEN;
SET <NUMBER> 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;
/* identify if it is a partial year wrt to IMS data - check if there are 12 months of data*/
num month_count=card(MONTH);
print month_count;
number launch_date{regimen};
read data regimen_Launch_dates into [regimen] launch_date;
print launch_date;
/*DEFINE THE PARAMETERS*/
var ps {REGIMEN, MONTH} init 0.3 >=0.00 <=.900,
annualpatients init 10000 <=25000 >=1000,
compliance {drug} >=0.3 <=.95,
dd{drug} >=4 <=6;
var factor >=0.2 <=1;
for {d in drug} fix dd[d]=dailydose[d];
for {d in drug} unfix dd["Ribavarin"];
if month_count=12 then fix factor=1;
else unfix factor;
impvar su_add{ d in drug, M IN MONTH}=units[d,m]/dd[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]/month_count;
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{r in regimen, m in month} (ps[r,M]*annualpatients/(avg_dot/(30/7)));
impvar total_units=sum{m in month,d in drug} su_add[d,M];
impvar monthly_units{m in month}=sum{d in drug} su_add[d,M];
impvar ratio{m in month}=monthly_units[m]/total_units;
impvar monthlypatients{m in month}=pat_tot*ratio[m]*factor;
/* 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]*dd[d]*compliance[d]*30*drug_share[d,m];
/* deviation*/
impvar deviation{d in drug, M IN MONTH}=(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}: -0.05 <= deviation[d,m] <= 0.05;
/*SET THE OBJECTIVE*/
var error_m_d {d in drug,M IN MONTH};
CON error_con {d in drug,M IN MONTH}: error_m_d[d,m] = (units[D,M]-calc_units[d,m])/dd[d];
Min error = sum{d in drug,M IN MONTH}error_m_d[D,M]**2;
/*min error{d in drug, m in month} = deviation[d,m];*/
/* expand;*/
solve with nlp;
create data pat_solver_shares from [regimen month] ps;
create data matching from [drug month] units calc_units deviation error_m_d;
create data compliance from [drug] compliance DD;
create data dot from [regimen]dot;
create data annual_pats from annualpatients;
create data factor from factor;
quit;
```

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to dpullarkat

02-23-2016 11:50 AM

Please also include the input data.

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to RobPratt

02-24-2016 12:26 AM

Hi Robb,

Below is the code, and attached are the input data file and also the output data from the optmodel procedure.

Thanks,

Deepthi

```
%let folder=C:\Users\pullarkd\Documents\Work\Solvers SAS Trial\;
%let input_file=Masked Data.xlsx;
libname solver "&folder.";
/*read in data from excel input template*/
LIBNAME inputs EXCEL PATH="&folder.\&input_file.";
data units_data_raw;
set inputs.'Inputs-Units Data$'n;
run;
data drug_consumption_master;
set inputs.drugdosing;
where dailydose is not null;
run;
data regimen_dot_master;
set inputs.regimen_dot_master;
where dot is not null;
run;
data regimen_drug_master;
set inputs.mapping;
where weight is not null;
run;
data compliance_constraints;
set inputs.compliance_range;
where upper is not null;
run;
data regimen_key_drug;
set inputs.regimen_key_drug;
where drug is not null;
run;
LIBNAME inputs CLEAR;
proc sql;
create table units_data_3 as
select country, period as month, drug, sum(value) as units
from units_data_raw
group by 1,2,3;
run;
proc sql;
create table month as
select distinct month
from units_data_3;
quit;
proc sql;
create table master as
select drug, month
from (select distinct drug from units_data_3) as a, month;
quit;
proc sql;
create table units_data as
select b.country, a.drug, a.month, case when units is null then 0 else units end as units
from master as a
left join units_data_3 as b
on a.drug=b.drug and a.month=b.month;
quit;
proc sql;
create table drug_consumption as
select * from drug_consumption_master
where drug in (select distinct drug from units_data);
quit;
proc sql;
create table regimen_drug as
select regimen, drug, weight
from regimen_drug_master
where regimen in
(select distinct regimen from
regimen_key_drug
where drug in (select distinct drug from units_data));
quit;
proc sql;
create table regimen_dot as
select distinct regimen, dot
from regimen_dot_master
where regimen in
(select distinct regimen from
regimen_key_drug
where drug in (select distinct drug from units_data));
quit;
proc sql;
create table launch_dates as
select drug, min(month) as launch_date
from
(select drug, month, units
from units_data where units >0) as a
group by 1;
quit;
proc sql;
create table regimen_launch_dates as
select regimen, min(launch_date) as launch_date
from
(select b.regimen, b.drug, a.launch_date, b.weight
from launch_dates as a
left join regimen_drug as b
on a.drug=b.drug) as a
where weight>0 and
(drug in ("Drug3", "Drug4", "Drug5", "Drug6", "Drug7", "Drug8", "Drug9") or regimen="Reg1")
group by 1;
quit;
PROC OPTMODEL;
SET <STRING> DRUG;
SET <STRING> REGIMEN;
SET <NUMBER> 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;
/*identify if it is a partial year wrt to IMS data - check if there are 12 months of data*/
num month_count=card(MONTH);
print month_count;
number launch_date{regimen};
read data regimen_Launch_dates into [regimen] launch_date;
/*DEFINE THE PARAMETERS*/
var ps {REGIMEN, MONTH} init 0.3 >=0.00 <=.900,
annualpatients init 10000 <=25000 >=1000,
dd{drug} >=4 <=6;
var factor >=0.2 <=1;
var compliance{drug};
read data compliance_constraints into [drug] compliance.lb=lower compliance.ub=upper;
for {d in drug} fix dd[d]=dailydose[d];
for {d in drug} unfix dd["Drug2"];
FOR {R IN REGIMEN,M IN MONTH} IF M<LAUNCH_DATE[r] THEN FIX PS[R,M]=0;
if month_count=12 then fix factor=1;
else unfix factor;
impvar su_add{ d in drug, M IN MONTH}=units[d,m]/dd[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]/month_count;
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{r in regimen, m in month} (ps[r,M]*annualpatients/(avg_dot/(30/7)));
impvar total_units=sum{m in month,d in drug} su_add[d,M];
impvar monthly_units{m in month}=sum{d in drug} su_add[d,M];
impvar ratio{m in month}=monthly_units[m]/total_units;
impvar monthlypatients{m in month}=pat_tot*ratio[m]*factor;
/* 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]*dd[d]*compliance[d]*30*drug_share[d,m];
/* deviation*/
impvar deviation{d in drug, M IN MONTH}=(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}: -0.05 <= deviation[d,m] <= 0.05;
/*SET THE OBJECTIVE*/
var error_m_d {d in drug,M IN MONTH};
CON error_con {d in drug,M IN MONTH}: error_m_d[d,m] = (units[D,M]-calc_units[d,m])/dd[d];
Min error = sum{d in drug,M IN MONTH}error_m_d[D,M]**2;
/*min error{d in drug, m in month} = deviation[d,m];*/
/* expand;*/
solve with nlp/maxiter=100;
create data pat_solver_shares from [regimen month] ps;
create data matching from [drug month] units calc_units deviation error_m_d;
create data compliance from [drug] compliance DD;
create data dot from [regimen]dot;
create data annual_pats from annualpatients;
create data factor from factor;
quit;
```

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to dpullarkat

02-25-2016 04:33 PM

You have stopped the solver early by setting MAXITER=100, and the resulting Infeasibility value is positive, meaning that some bounds or constraints are violated. Are you sure that the problem has a feasible solution?

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to RobPratt

02-26-2016 12:13 AM

Hi Robb,

Sorry, I had forgotten to remove the MAXITER option from the code I shared (it was for testing purposes).

But the issue is still present even if the optmodel ran 5000 iterations, which is the maximum limit (or can I increase it using the MAXITER?).

The problem does not have a 'fully' feasible solution in all scenarios, but it is more important that the variables are within the constraints, even if the objective value is somewhat large. How does the proc optmodel algorithm work? Does is use all possible values of the variables, even those outside the constraints?

If I reduce the scale of the objective, by mutiplying with e.g. 10e-4, would it work?

Thanks,

Deepthi

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to dpullarkat

02-26-2016 04:34 PM

There are two main SAS nlp solvers available. One is IP (interior point solver), and the other is AS (active set solver). Currently, constraints, including variable bounds are relaxed, in solver IP, in order to permit larger steps to enhance the performance. When a feasible point is not found, some of the constraints can be violated, including bound constraints. On the other hand, solver AS will honor variable bounds. So if you prefer the algorithm stays inside the box bounds, for now, you might should use solver AS.

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to wezhou

02-29-2016 07:57 AM

Thanks for helping out!

Is there any good material that can help me in understanding the mathematical basis for this? I was not able to find the information I wanted from support.sas

Thanks again,

Deepthi

- Mark as New
- Bookmark
- Subscribe
- RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to dpullarkat

02-29-2016 04:28 PM

I hope the following links can be helpful.

http://www.neos-guide.org/content/gradient-projection-methods

http://www.ccom.ucsd.edu/~peg/papers/pdmerit.pdf