RobPratt Tracker

Re: Proc Optmodel: How to find the top two or three solutions?

RobPratt — Mon, 12 Jan 2026 19:07:45 GMT

You can use the MILP solver options MAXPOOLSOLS=3 and SOLTYPE=best to find the best 3 solutions:

A less efficient but more flexible approach is to introduce dynamically generated no-good constraints to cut off each solution as it is found:

   num numSolutions init 0;
   set SOLUTIONS = 1..numSolutions;
   set <str,num,num> SUPPORT_sol {SOLUTIONS};
   con NoGood {sol in SOLUTIONS}:
      sum {<c,a,e> in SUPPORT_sol[sol]} (1 - CHI[c,a,e]) >= 1;

   for {1..3} do;
      solve;
      numSolutions = numSolutions + 1;
      SUPPORT_sol[numSolutions] = {<c,a,e> in stations cross orientations: CHI[c,a,e] > 0.5};
   end;

   create data work.OptimalSolution_RoO from [sol station a e]={sol in SOLUTIONS, <c,a,e> in SUPPORT_sol[sol]} CHI=1;

In the FOR loop, you could alternatively replace the SOLVE statement with your more intricate solution process that loops over pairs to find a good solution.

Re: Proc Optmodel - output

RobPratt — Fri, 28 Nov 2025 16:57:01 GMT

The closest you can get to returning a BestBound solution is to solve the LP relaxation:

solve relaxint;

The resulting objective value will often be close to the first or second BestBound value that appears in the iteration log. The first BestBound value sometimes arises from a trivial bound that considers only the objective function and the variable bounds.

Re: Proc Optmodel - output

RobPratt — Fri, 28 Nov 2025 16:31:03 GMT

What you have is correct, except that the RUN should instead be QUIT.

The values of H and X in the resulting solution data set will correspond to the objective value reported in the BestInteger column. You will see the same values if you execute this statement:

print H X;

You should also see the value in BestInteger if you execute this statement:

print obj;

What makes you suspect that the values of H and X instead correspond to the BestBound column?

Re: Proc Optmodel Result: Why is Objective = 0?

RobPratt — Tue, 14 Oct 2025 15:55:24 GMT

When there are two stations and each one can use at most one orientation, it is not possible for a target to be triply covered, so IsTriplyCovered[t] = 0 for all t, and the maximum objective value is 0. There are multiple optimal solutions, and the one returned by the solver has only one CHI variable > 0.5.

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?!

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 2 Find the BEST Improvement on DIF between _tmin_orig and _tmin_x[treatment/method] for BLOC

RobPratt — Thu, 09 Oct 2025 16:55:36 GMT

Please open a new question for that.

Re: How 2 Find the BEST Improvement on DIF between _tmin_orig and _tmin_x[treatment/method] for BLOC

RobPratt — Fri, 03 Oct 2025 21:08:55 GMT

TotalScore.sol is the value of TotalScore from the most recent solver call. You can alternatively use a manually specified threshold. The following code changes yield the same objective values as in your Excel file:

   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 PUT statement shows the resulting objective values:

TotalScore=-3 TotalDif=13365

Re: How 2 Find the BEST Improvement on DIF between _tmin_orig and _tmin_x[treatment/method] for BLOC

RobPratt — Thu, 02 Oct 2025 22:45:02 GMT

I see a few issues here:

You cannot have = expression in a VAR statement. You should keep the original impvar Score and max TotalScore declarations.
The outcome parameter requires four arguments, not two: outcome[d,b,t,m], not outcome[d,b].
You have not defined _time_orig[b]. Is this maybe supposed to be _tmin_orig[d,b,t]?
Your MinScore declaration is missing a colon, and the right-hand side should depend on the previous objective function value. Correct syntax is: con MinScore: TotalScore >= TotalScore.sol;

I'm not quite sure what you want the secondary objective to be, but here's an outline of the two stages:

   /* primary objective */
   solve;

   con MinScore: TotalScore >= TotalScore.sol;

   impvar difx {<d,b> in DATE_BLOCK} = sum {<(d),(b),t> in DATE_BLOCK_TREATMENT, m in METHODS} (_tmin_orig[d,b,t] - outcome[d,b,t,m]) * SelectTreatmentMethod[t,m];
   max TotalDif = sum {<d,b> in DATE_BLOCK} difx[d,b];

   /* secondary objective */
   solve with milp / primalin;

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Wed, 01 Oct 2025 14:43:45 GMT

I would be glad to help, but first please create a new question for this part.

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Mon, 29 Sep 2025 13:46:30 GMT

   con CardinalityTreatment234:
      sum {t in 2..4} SelectTreatment[t] = 1;

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Tue, 23 Sep 2025 04:52:53 GMT

Correct syntax is:

   create data want_x_dtbloc from [d b]=DATE_BLOCK;
   create data want_trtmt from [t]=TREATMENTS;

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Tue, 23 Sep 2025 15:40:15 GMT

To match the Excel logic, you can change the optimization model as follows:

   /* define optimization model */
   var SelectTreatment {TREATMENTS} binary;
   var SelectMethod {METHODS} binary;
   var SelectTreatmentMethod {TREATMENTS, METHODS} binary;
   var IsGood {DATE_BLOCK} binary;
   /* if IsGood[d,b] = 1 then Score[d,b] = 1 else Score[d,b] = -5 */
   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 CardinalityTreatment:
      sum {t in TREATMENTS} SelectTreatment[t] = 2;
   con CardinalityMethod:
      sum {m in METHODS} SelectMethod[m] = 1;
   /* SelectTreatmentMethod[t,m] <= SelectTreatment[t] * SelectMethod[m] */
   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];

The idea is that if IsGood[d,b] = 1, the AtLeastOneGood constraint will force SelectTreatmentMethod[t,m] = 1 for at least one t and m such that outcome[d,b,t,m] < _tmin_orig[d,b,t]. And the Linearize1 and Linearize2 constraints together force SelectTreatment[t] = 1 and SelectMethod[m] = 1 for that t and m.

With these changes, the optimal objective value is 3, and there are 13 such optimal solutions for SelectTreatment.

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Sat, 20 Sep 2025 20:42:14 GMT

Because the two cardinality constraints are equalities, increasing the right-hand side values does not make the problem looser. If the = were instead <=, then increasing the right-hand side would indeed yield a relaxation. Do you want to select exactly the specified numbers or at most the specified numbers?

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Thu, 18 Sep 2025 14:53:37 GMT

It looks like you are using an old version of SAS. That dynamic data set name functionality was introduced in SAS/OR 14.1 in SAS 9.4M3 in 2015:

Data Set Input/Output :: SAS/OR(R) 14.1 User's Guide: Mathematical Programming

Here's a workaround that uses PROC APPEND inside a SUBMIT block:

   /* call MILP solver to find top few solutions for SelectTreatment */
   num numSolsWanted = 5;
   num numSolsFound init 0;
   set SOLS = 1..numSolsFound;
   set TREATMENTS_s {SOLS};
   con ExcludeSolution {s in SOLS}:
      sum {t in TREATMENTS_s[s]} SelectTreatment[t] <= card(TREATMENTS_s[s]) - 1;
   for {1..numSolsWanted} do;
      put numSolsFound=;
      solve;
      numSolsFound = numSolsFound + 1;
      TREATMENTS_s[numSolsFound] = {t in TREATMENTS: SelectTreatment[t].sol > 0.5};
      create data want_method from [method] sol=numSolsFound SelectMethod;
      create data want_treatment from [treatment] sol=numSolsFound SelectTreatment;
      create data want_dt_bloc from [dt bloc] sol=numSolsFound IsGood Score;
      submit;
         proc append base=want_method_all data=want_method;
         run;
         proc append base=want_treatment_all data=want_treatment;
         run;
         proc append base=want_dt_bloc_all data=want_dt_bloc;
         run;
      endsubmit;
   end;

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Thu, 18 Sep 2025 00:36:34 GMT

To get the 5 best pairs of treatments, you can replace the SOLVE and CREATE DATA statements as follows:

   /* call MILP solver to find top few solutions for SelectTreatment */
   num numSolsWanted = 5;
   num numSolsFound init 0;
   set SOLS = 1..numSolsFound;
   set TREATMENTS_s {SOLS};
   con ExcludeSolution {s in SOLS}:
      sum {t in TREATMENTS_s[s]} SelectTreatment[t] <= card(TREATMENTS_s[s]) - 1;
   for {1..numSolsWanted} do;
      put numSolsFound=;
      solve;
      numSolsFound = numSolsFound + 1;
      TREATMENTS_s[numSolsFound] = {t in TREATMENTS: SelectTreatment[t].sol > 0.5};
      create data ('want_method'||numSolsFound) from [method] SelectMethod;
      create data ('want_treatment'||numSolsFound) from [treatment] SelectTreatment;
      create data ('want_dt_bloc'||numSolsFound) from [dt bloc] IsGood Score;
   end;

Re: HOW TO Find the Best Combination Of Treatments for Date/BLOC[dt/bloc] ?!

RobPratt — Wed, 17 Sep 2025 19:52:42 GMT

If I understand correctly, the following code does what you want:

data have;
   set mm_out_x6_x;
   rename _tmin=method1 _tmin_ew=method2 _tmin_wt_r=method3 _tmin_wt_r2=method4;
run;

proc optmodel;
   /* read input data */
   set METHODS = 1..4;
   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 SelectMethod {METHODS} binary;
   var SelectTreatment {TREATMENTS} binary;
   var IsGood {DATE_BLOCK} binary;
   /* if IsGood[d,b] = 1 then Score[d,b] = 1 else Score[d,b] = -5 */
   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 CardinalityMethod:
      sum {m in METHODS} SelectMethod[m] = 1;
   con CardinalityTreatment:
      sum {t in TREATMENTS} SelectTreatment[t] = 2;
   /* if IsGood[d,b] = SelectMethod[m] = SelectTreatment[t] = 1 then outcome[d,b,t,m] < _tmin_orig[d,b,t] */
   con NoGood {<d,b,t> in DATE_BLOCK_TREATMENT, m in METHODS: outcome[d,b,t,m] >= _tmin_orig[d,b,t]}:
      IsGood[d,b] + SelectMethod[m] + SelectTreatment[t] <= 2; 

   /* call MILP solver */
   solve;

   /* create output data */
   create data want_method from [method] SelectMethod;
   create data want_treatment from [treatment] SelectTreatment;
   create data want_dt_bloc from [dt bloc] IsGood Score;
quit;

An optimal solution selects method 3 and treatments 12 and 15, yielding a maximum total score of -21.

If this does not match your expectation, please provide a sample solution, together with the calculation of its total score.

Re: How to find the BEST combination of two CID to Score this?!

RobPratt — Tue, 16 Sep 2025 19:40:38 GMT

Here's a way to use the runOptmodel action with groupBy to solve a separate problem for each dt and bloc:

data have;
   set lib.mmind_temp_x;
   if _tmin_orig > _tmin OR _tmin_orig > _tmin_ew OR _tmin_orig > _tmin_wt_r OR _tmin_orig > _tmin_wt_r2
   then score = 1;
   else score = -5;
run;

data mycas.have;
   set have;
run;

proc cas;
   source pgm;
      /* read input data */
      set IDS;
      num score {IDS};
      read data have into IDS=[condi_id] score;

      /* define optimization model */
      var Select {IDS} binary;
      max TotalScore = sum {i in IDS} score[i] * Select[i];
      con Cardinality:
         sum {i in IDS} Select[i] = 2;

      /* call MILP solver */
      solve;

      /* create output data */
      create data want from [condi_id] score Select;
   endsource;

   action optimization.runOptmodel / code=pgm groupBy={'dt','bloc'};
quit;

But you can get the same results by sorting in descending order of score and keeping the top two:

proc sort data=have out=want;
   by dt bloc descending score;
run;

data want;
   set want;
   by dt bloc;
   retain rank;
   if first.bloc then rank = 1;
   else rank + 1;
   Select = (rank <= 2);
run;

proc sort data=want;
   by dt bloc condi_id;
run;

Re: How to find the BEST combination of two CID to Score this?!

RobPratt — Tue, 16 Sep 2025 15:55:18 GMT

To help clarify the problem, can you please share a feasible (not necessarily optimal) solution and the corresponding objective values?

Re: find all possible sum of variables that amount to a value

RobPratt — Wed, 10 Sep 2025 14:49:57 GMT

@Ksharp Sorry, I had oversimplfiied the code, which I have now corrected.

The corrected line is:

set VARS {IDS} init VARS_ALL;