topic Re: How to build/code up a penalty function?! in Mathematical Optimization, Discrete-Event Simulation, and OR

How to build/code up a penalty function?!

hellohere — Fri, 10 Oct 2025 02:01:19 GMT

@RobPratt

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=;

Re: How to build/code up a penalty function?!

hellohere — Fri, 10 Oct 2025 02:02:18 GMT

Initial Thread and dataset: https://communities.sas.com/t5/Mathematical-Optimization/How-2-Find-the-BEST-Improvement-on-DIF-between-tmin-orig-and/m-p/976238#M4381

Re: How to build/code up a penalty function?!

RobPratt — Fri, 10 Oct 2025 14:19:51 GMT

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];

Re: How to build/code up a penalty function?!

hellohere — Sat, 11 Oct 2025 02:36:15 GMT

@RobPratt

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;

Re: How to build/code up a penalty function?!

hellohere — Sat, 11 Oct 2025 02:54:00 GMT

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.

Re: How to build/code up a penalty function?!

hellohere — Sat, 11 Oct 2025 03:26:35 GMT

@RobPratt

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].

Re: How to build/code up a penalty function?!

RobPratt — Tue, 14 Oct 2025 15:03:01 GMT

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);

Re: How to build/code up a penalty function?!

hellohere — Tue, 14 Oct 2025 00:03:19 GMT

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.