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 (but with -20 instead of -10) 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);
