How to code-up a penalty function on the difference between _tmin_orig and _tmin_x?!
If _tmin_orig<=_tmin_x, score=-10
else if _tmin_orig>_tmin_x+60, score=2
else if _tmin_orig>_tmin_x+36, score=4
else if _tmin_orig>_tmin_x+24, score=8
else if _tmin_orig>_tmin_x+12, score=6
else if _tmin_orig>_tmin_x+6, score=2
else if _tmin_orig>_tmin_x+0, score=1
Thanks,
con CardinalityMethod:
sum {m in METHODS} SelectMethod[m] = 1;
/* primary objective */
solve;
/* con MinScore: TotalScore >= TotalScore.sol;*/
con MinScore: TotalScore >= -5;
impvar difx {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m];
max TotalDif = sum {<d,b> in DATE_BLOCK} difx[d,b];
/* secondary objective */
solve with milp / primalin;
put TotalScore= TotalDif=;
The Excel file makes it clearer what you are trying to do.
Here's a way to use the black-box solver:
impvar Difx {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m];
impvar Penalty {<d,b> in DATE_BLOCK} =
if DifX[d,b] <= 0 then -20
else if DifX[d,b] > 60 then 2
else if DifX[d,b] > 36 then 4
else if DifX[d,b] > 24 then 8
else if DifX[d,b] > 12 then 6
else if DifX[d,b] > 6 then 2
else if DifX[d,b] > 0 then 1;
/* impvar Penalty {<d,b> in DATE_BLOCK} = penaltyCoefficient(DifX[d,b]);*/
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
/* secondary objective */
solve obj TotalPenalty with blackbox / primalin;
Note that the commented out line instead uses the FCMP function from my earlier reply and is mathematically equivalent.
But the solution returned by the black-box solver is not necessarily globally optimal. Here's an alternative approach that uses the MILP solver to find a globally optimally solution:
num upperBoundDifx =
max {<d,b> in DATE_BLOCK}
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0);
num epsilon = 1e-2;
num numIntervals = 7;
set INTERVALS = 1..numIntervals;
num ub {INTERVALS};
ub[1] = 0;
ub[2] = 6;
ub[3] = 12;
ub[4] = 24;
ub[5] = 36;
ub[6] = 60;
ub[7] = upperBoundDifx;
num lb {i in INTERVALS} = (if i = 1 then 0 else ub[i-1] + epsilon);
num p {INTERVALS} = [-20, 1, 2, 6, 8, 4, 2];
print lb ub p;
var IsInterval {DATE_BLOCK, INTERVALS} binary;
con OneInterval {<d,b> in DATE_BLOCK}:
sum {i in INTERVALS} IsInterval[d,b,i] = 1;
con LeftCon {<d,b> in DATE_BLOCK}:
sum {i in INTERVALS} lb[i] * IsInterval[d,b,i] <= DifX[d,b];
con RightCon {<d,b> in DATE_BLOCK}:
sum {i in INTERVALS} ub[i] * IsInterval[d,b,i] >= DifX[d,b];
impvar Penalty {<d,b> in DATE_BLOCK} = sum {i in INTERVALS} p[i] * IsInterval[d,b,i];
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
/* secondary objective */
solve with milp / primalin;
If you know that the values of DifX will always take integer values, you can use 1 instead of 1e-2 for epsilon.
For your input data, both black-box and MILP return a solution with TotalPenalty = 2.
If you fix the selected treatments and methods to the ones selected in the Excel file and then call either solver, the resulting solutions match the Excel file and yield TotalPenalty = -7:
/* temporarily fix to match Excel solution */
for {t in TREATMENTS} fix SelectTreatment[t] = (t in {13,15});
for {m in METHODS} fix SelectMethod[m] = (m = 3);
Assuming you want to minimize the total penalty, here is one way:
impvar Penalty {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
if _tmin_orig[d,b,t] - outcome[d,b,t,m] <= 0 then -10
else if _tmin_orig[d,b,t] - outcome[d,b,t,m] > 60 then 2
else if _tmin_orig[d,b,t] - outcome[d,b,t,m] > 36 then 4
else if _tmin_orig[d,b,t] - outcome[d,b,t,m] > 24 then 8
else if _tmin_orig[d,b,t] - outcome[d,b,t,m] > 12 then 6
else if _tmin_orig[d,b,t] - outcome[d,b,t,m] > 6 then 2
else if _tmin_orig[d,b,t] - outcome[d,b,t,m] > 0 then 1
) * SelectTreatmentMethod[t,m];
min TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
An alternative approach is to use an FCMP function:
proc fcmp outlib=work.funcs.test;
function penaltyCoefficient(dif);
if dif <= 0 then return(-10);
else if dif > 60 then return(2);
else if dif > 36 then return(4);
else if dif > 24 then return(8);
else if dif > 12 then return(6);
else if dif > 6 then return(2);
else if dif > 0 then return(1);
endsub;
run;
option cmplib=work.funcs;
data try;
do x = -5 to 70 by 5;
p = penaltyCoefficient(x);
output;
end;
run;
Then the Penalty declaration can be written more compactly:
impvar Penalty {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} penaltyCoefficient(_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m];
Thanks for helps.
I tried the codes below, two variations. Both complains "ERROR: The specified optimization technique does not allow nonlinear objectives." .
Note, TotalPenalty is like TotalScore. It is to MAX'it. Also -10 is swtiched to -20.
With Treatment/Method, if _tmin_x is less than _tmin_orig, improvement[Score>0].
The Score is consider by each BLOC/DT.
impvar Penalty {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) <= 0 then -20
else if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) > 60 then 2
else if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) > 36 then 4
else if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) > 24 then 8
else if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) > 12 then 6
else if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) > 6 then 2
else if sum((_tmin_orig[d,b,t] - outcome[d,b,t,m])* SelectTreatmentMethod[t,m]) > 0 then 1
) * SelectTreatmentMethod[t,m];
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
impvar Penalty {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] <= 0 then -20
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 60 then 2
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 36 then 4
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 24 then 8
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 12 then 6
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 6 then 2
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 0 then 1
) * SelectTreatmentMethod[t,m];
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
Complete codes are below. Also attached is the dataset.
%let optds=Mm_out_x6_x; %let methct=4;
data have;
set &optds.;
rename _tmin=method1 _tmin_ew=method2 _tmin_wt_r=method3 _tmin_wt_r2=method4;
run;
title "OPT on &optds.";
proc optmodel;
/* read input data */
set METHODS = 1..&methct.;
set <str,num,num> DATE_BLOCK_TREATMENT;
num _tmin_orig {DATE_BLOCK_TREATMENT};
num outcome {DATE_BLOCK_TREATMENT, METHODS};
read data have into DATE_BLOCK_TREATMENT=[dt bloc condi_id] _tmin_orig
{m in METHODS} <outcome[dt,bloc,condi_id,m]=col('method'||m)>;
set TREATMENTS = setof {<d,b,t> in DATE_BLOCK_TREATMENT} t;
set DATE_BLOCK = setof {<d,b,t> in DATE_BLOCK_TREATMENT} <d,b>;
/* define optimization model */
var SelectTreatment {TREATMENTS} binary;
var SelectMethod {METHODS} binary;
var SelectTreatmentMethod {TREATMENTS, METHODS} binary;
var IsGood {DATE_BLOCK} binary;
/* IT WAS*/
impvar Score {<d,b> in DATE_BLOCK} = 1 * IsGood[d,b] - 5 * (1 - IsGood[d,b]);
max TotalScore = sum {<d,b> in DATE_BLOCK} Score[d,b];
/*con CardinalityTrtmtBloc:
sum {t in 5..6} SelectTreatment[t] = 1; /*enforce to select ONE AND ONLY ONE among 2-4[or alike]*/
con CardinalityTreatment:
sum {t in TREATMENTS} SelectTreatment[t] = 2;
con CardinalityMethod:
sum {m in METHODS} SelectMethod[m] = 2;
con Linearize1 {t in TREATMENTS, m in METHODS}:
SelectTreatmentMethod[t,m] <= SelectTreatment[t];
con Linearize2 {t in TREATMENTS, m in METHODS}:
SelectTreatmentMethod[t,m] <= SelectMethod[m];
con AtLeastOneGood {<d,b> in DATE_BLOCK}:
IsGood[d,b] <= sum {<(d),(b),t> in DATE_BLOCK_TREATMENT,
m in METHODS: outcome[d,b,t,m] < _tmin_orig[d,b,t]} SelectTreatmentMethod[t,m];
/* call MILP solver , TOP one only
solve;
create data want_method from [method] SelectMethod;
create data want_treatment from [treatment] SelectTreatment;
create data want_dt_bloc from [dt bloc] IsGood Score;
create data want_x_dtbloc from [d b]=DATE_BLOCK;
create data want_trtmt from [t]=TREATMENTS;*/
/* source: https://communities.sas.com/t5/Mathematical-Optimization/How-2-Find-the-BEST-Improvement-on-DIF-between-tmin-orig-and/m-p/976283/highlight/false#M4384
*/
/* primary objective */
solve;
create data want_method from [method] SelectMethod;
create data want_treatment from [treatment] SelectTreatment;
create data want_dt_bloc from [dt bloc] IsGood Score;
create data want_x_dtbloc from [d b]=DATE_BLOCK;
create data want_trtmt from [t]=TREATMENTS;
/* con MinScore: TotalScore >= TotalScore.sol;*/
con MinScore: TotalScore >= -5;
/*impvar difx {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m];
max TotalDif = sum {<d,b> in DATE_BLOCK} difx[d,b];
*/
impvar Penalty {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] <= 0 then -20
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 60 then 2
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 36 then 4
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 24 then 8
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 12 then 6
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 6 then 2
else if max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m] > 0 then 1
) * SelectTreatmentMethod[t,m];
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
/* secondary objective */
solve with milp / primalin;
put TotalScore= TotalPenalty=;
create data want_method_f from [method] SelectMethod;
create data want_treatment_f from [treatment] SelectTreatment;
create data want_dt_bloc_f from [dt bloc] IsGood Score;
create data want_x_dtbloc_f from [d b]=DATE_BLOCK;
create data want_trtmt_f from [t]=TREATMENTS;
quit;
I also tried the code below, which I bet is more meaningful.
impvar Penalty {<d,b> in DATE_BLOCK} =
if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] < 0 then -20
else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 60 then 2
else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 36 then 4
else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 24 then 8
else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 12 then 6
else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 6 then 2
else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] >= 0 then 1;
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
It still has the same complain. Might need to set nonlinear approach?!
21702 impvar Penalty {<d,b> in DATE_BLOCK} =
21703 if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21704 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] < 0 then -20
21705 else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21706 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 60 then 2
21707 else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21708 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 36 then 4
21709 else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21710 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 24 then 8
21711 else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21712 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 12 then 6
21713 else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21714 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] > 6 then 2
21715 else if sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (
21716 _tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m] >= 0 then 1;
21717
21718 max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
21719
21720 /* secondary objective */
21721 solve with milp / primalin;
NOTE: Problem generation will use 4 threads.
NOTE: The problem has 109 variables (0 free, 0 fixed).
NOTE: The problem uses 31 implicit variables.
NOTE: The problem has 109 binary and 0 integer variables.
NOTE: The problem has 162 linear constraints (159 LE, 2 EQ, 1 GE, 0 range).
NOTE: The problem has 819 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
ERROR: The specified optimization technique does not allow nonlinear objectives.
Attached is the EXCEL/SOLVER version. The Penalty is only the function on difx[from previous].
Note the columns AE:AJ are completely identical to the columns U:Z, because The Penalty is only the function on difx[from previous].
The Excel file makes it clearer what you are trying to do.
Here's a way to use the black-box solver:
impvar Difx {<d,b> in DATE_BLOCK} =
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0) * SelectTreatmentMethod[t,m];
impvar Penalty {<d,b> in DATE_BLOCK} =
if DifX[d,b] <= 0 then -20
else if DifX[d,b] > 60 then 2
else if DifX[d,b] > 36 then 4
else if DifX[d,b] > 24 then 8
else if DifX[d,b] > 12 then 6
else if DifX[d,b] > 6 then 2
else if DifX[d,b] > 0 then 1;
/* impvar Penalty {<d,b> in DATE_BLOCK} = penaltyCoefficient(DifX[d,b]);*/
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
/* secondary objective */
solve obj TotalPenalty with blackbox / primalin;
Note that the commented out line instead uses the FCMP function from my earlier reply and is mathematically equivalent.
But the solution returned by the black-box solver is not necessarily globally optimal. Here's an alternative approach that uses the MILP solver to find a globally optimally solution:
num upperBoundDifx =
max {<d,b> in DATE_BLOCK}
sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} max(_tmin_orig[d,b,t] - outcome[d,b,t,m], 0);
num epsilon = 1e-2;
num numIntervals = 7;
set INTERVALS = 1..numIntervals;
num ub {INTERVALS};
ub[1] = 0;
ub[2] = 6;
ub[3] = 12;
ub[4] = 24;
ub[5] = 36;
ub[6] = 60;
ub[7] = upperBoundDifx;
num lb {i in INTERVALS} = (if i = 1 then 0 else ub[i-1] + epsilon);
num p {INTERVALS} = [-20, 1, 2, 6, 8, 4, 2];
print lb ub p;
var IsInterval {DATE_BLOCK, INTERVALS} binary;
con OneInterval {<d,b> in DATE_BLOCK}:
sum {i in INTERVALS} IsInterval[d,b,i] = 1;
con LeftCon {<d,b> in DATE_BLOCK}:
sum {i in INTERVALS} lb[i] * IsInterval[d,b,i] <= DifX[d,b];
con RightCon {<d,b> in DATE_BLOCK}:
sum {i in INTERVALS} ub[i] * IsInterval[d,b,i] >= DifX[d,b];
impvar Penalty {<d,b> in DATE_BLOCK} = sum {i in INTERVALS} p[i] * IsInterval[d,b,i];
max TotalPenalty = sum {<d,b> in DATE_BLOCK} Penalty[d,b];
/* secondary objective */
solve with milp / primalin;
If you know that the values of DifX will always take integer values, you can use 1 instead of 1e-2 for epsilon.
For your input data, both black-box and MILP return a solution with TotalPenalty = 2.
If you fix the selected treatments and methods to the ones selected in the Excel file and then call either solver, the resulting solutions match the Excel file and yield TotalPenalty = -7:
/* temporarily fix to match Excel solution */
for {t in TREATMENTS} fix SelectTreatment[t] = (t in {13,15});
for {m in METHODS} fix SelectMethod[m] = (m = 3);
Thanks a lot. I bet you take great efforts to help. I guess I need to pick up own IML.
No need to be global optimal.
It's finally time to hack! Remember to visit the SAS Hacker's Hub regularly for news and updates.