Min-max objective and max operator with concave function - AMPL with gurobi

[Syed Hasan]

Dear all,

I am running my model written in AMPL using MILP solver Gurobi at the
NEOS site.

Background:
----------------------------------------------------------------
Original cost function in the model: linear:

var C {i in NODES} = I[i] * c_f + (c_b * alpha * V[i]);

Objective function:
minimize cost: sum {i in NODES} C[i];

where 'I[i]' is indicator variable (binary) (0: node not selected, 1:
node selected)
'V[i]' is a variable defining traffic delivered from that node to all
other nodes. (I do have an expression for it)
c_f, c_b are fixed values i.e., 10 and alpha = 0.7


Now, I want to consider two cost variants:
---------------------------------------------------------------------

I. Non-linear cost component in the cost function:

Concave increasing function: (c_b * V(i))^g)
where g is 0.8.

Thus, I have incorporated the non-linear cost function in the AMPL
model by means of piece-wise linear approximation following AMPL's
notation (breakpoints, slopes) i.e., concave piecewise separable
function. 4-piece approximation with 3 bkpts and 4 slopes. Also,
formulating as min-max objective function due to concave piece-wise
linear terms.

var Z;
minimize max_cost: Z;
subject to cost_def:
Z >= sum {i in NODES} <<{p in 1..npce[i]-1} bkp[i,p];
{p in 1..npce[i]} slp[i,p]>> V[i] + sum {i in NODES} c_f *
I[i];


DATA:

bkpts:
1 1 10
1 2 20
1 3 30
.....
same bkpts sequence (10,20,30) for all 10 nodes.

slopes:
1 1 3.162
1 2 2.156
1 3 1.89
1 4 1.73

-- same slopes sequence for all 10 nodes

Slopes are computed from the function values and bkpts:
At V = 10, function value f(10) = (c_b*10)^0.75 = (10*10)^0.75 =
31.622

Thus slope between 0 (default) and 10: f(10)-f(0) / (bkpt1 - 0)
(31.622-0)/(10-0) = 3.1622 (piece-1)
similarly, for 10-20, 20-30..
and the last slope for the 4th piece is computed as: (f(40)-f(30)) /
(40 -bkpt3).
(I am not sure whether this is the right way).

ACTUAL PROBLEM/ISSUE:
----------------------------------------------------------------------------
The solver returned the result for the input model based on
constraints:
Only 1-node selected,i.e,'3', V[3] = 810, and the max_cost = 1431.48.
But I expect the answer to be following the concave function:
(810*10)^0.75 + 10 = 863.81

Other attempts:
-- I am considering that AMPL's delta-form is applied to my model for
the transformations.
-- Even I tried to add the lamda-formulation (inverted V-form
translation: each P-L term replaced by a convex combination of the
function values at the breakpoints (plus SOS-type 2 constraint)
---------------------------------------------------------------------------------------------
var w {i in NODES} = sum {p in 1..npce[i]-1} lamda[i,p] * (c_b *
bkp[i,p])^(0.75);

var x {i in NODES} = sum {p in 1..npce[i]-1} lamda[i,p] * bkp[i,p];

subject to lam {i in NODES}: sum {p in 1..npce[i]-1} lamda[i,p] = 1;

Adjacency constraint to drive SOS-type 2 i.e. either lamda[i,1] = 0 or
lamda[i,3] = 0; I haven't tried it in the model or at-least not done
correctly.

II. max (c_f, (c_b * V[i])^g)
---------------------------------------------------------------------------------
How to incorporate max operator with two args: parameter and concave
piecewise separable function in the objective funtion?. I have
followed from AMPL book 1993 edition (chapter 14: Piecewise linear
programs), but didn't got the way to solve for my problem, but I think
that some piece-wise linear terminology will have the answer for that
special piecewise linear function as written in II. AMPL FAQ does tell
about non-linearity due to max operator.

If you could provide me any guidance/help in my problems, I will be
very grateful.

Thanks!

Best Regards,
Syed

[Robert Fourer]

In variant I, it's not clear why you are defining the variable Z. You could
write

minimize max_cost:

sum {i in NODES} <<{p in 1..npce[i]-1} bkp[i,p];
{p in 1..npce[i]} slp[i,p]>> V[i]
+ sum {i in NODES} c_f * I[i];

In evaluating the result, keep in mind that AMPL assumes by default that a
piecewise-linear function takes the value zero when the variable is zero.
If you want it to be some other value at zero, add a constant to the
objective; if you want it to be zero at some other value, replace V[i] by
(V[i],...) where ... is an expression for the value at which you want the
function to be zero. These changes will not affect the optimal solution but
they do make a difference to the reported optimal value of the objective.

It is hard for me to offer more advice than this, as I do not follow your
example of the issue you are having. If V[3] = 810 but you only have
breakpoints at 10, 20, and 30, then you would not expect to get exactly the
same objective value from your piecewise-linear approximation as the value
from the actual nonlinear objective. The values should be exactly the same
only at the breakpoints.

In variant II, to minimize the maximum of two expressions, then it makes
sense to minimize Z, subject to Z being >= each of them.

Bob Fourer
4...@ampl.com

[Syed Hasan]

Hi Bob,

Thanks for your reply..
Cost variant I.
-- I had defined a new variable Z since I am approximating concave
function for which a maximization is needed i.e., a min-max objective.
But, with my realistic bkpts, I am getting the same result for both
schemes (Z (new variable) and minimize objective (default).
-- Now I am considering 7 pieces with 6 breakpts (135, 270,.....810)
and 7 slopes. The value 810 is the total traffic to be delivered. The
approximated objective value is close to the actual non-linear cost
function (objective), and with a difference (between 3-10).

Refs: Expressing Special Structures in an Algebraic Modeling Language
for Mathematical Programming and AMPL book (Chp:PWL programs).

Cost variant-II.
-- I am having difficulty in writing the max function with
expressions: c_f and piecewise linear concave function i.e., max (c_f,
(c_b * V[i])^g) in an AMPL model/Gurobi.
(a) I have followed this approach, but this might work only if c_f =
0, but in my case c_f is a parameter with fixed values: 10, 40, 100,
etc.

var Z {i in NODES} >= c_f;
minimize max_cost: sum {i in NODES} Z[i];
subject to Zref {i in NODES}:
Z[i] >= <<{p in 1..npce[i]-1} bkp[i,p];
{p in 1..npce[i]} slp[i,p]>> V[i];

(b) Also, I followed AMPL book (chp: 17, other PWL functions which
gave an example of max function with PWL penalty function), but didn't
quite able to re-write the max function for my case. Could you provide
your insights(re-writing, linearization) towards this problem ?.

Thanks!

- Syed

[Robert Fourer]

This looks like a proper linearization to me:

var Z {i in NODES} >= c_f;
minimize max_cost: sum {i in NODES} Z[i];
subject to Zref {i in NODES}:
Z[i] >= <<{p in 1..npce[i]-1} bkp[i,p];
{p in 1..npce[i]} slp[i,p]>> V[i];

Z[i] is >= c_f and also >= <<...;...>> V[i], so it must be >= the max of
those two. But I don't understand what you mean when you say, "c_f is a
parameter with fixed values: 10, 40, 100, etc." since a parameter can have
only one value at a time. If you want to say that Z[i] may only equal 10 or
40 or 100 or ..., then that is another matter, but then you would not have a
min-max problem.

--
You received this message because you are subscribed to the Google Groups
"AMPL Modeling Language" group.
To post to this group, send email to am...@googlegroups.com.
To unsubscribe from this group, send email to
ampl+uns...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/ampl?hl=en.

[Syed Hasan]

Hi Bob,

Thanks. regarding c_f, I had not written properly in my post, but I
meant to say c_f (unique value) for different runs (experiments) to
obtain different objective values. I am using the same linearization
as posted in my model, I will get back to you with the issues I am
facing with an example.

Thanks!
-- Syed

[Syed Hasan]

Bob,
In my case, the objective function also has the indicator variable
I[i], thus I do have some problem.

var Z {i in NODES} >= c_f;

minimize max_cost: sum {i in NODES} Z[i] * I[i];

subject to Zref {i in NODES}:
Z[i] >= <<{p in 1..npce[i]-1} bkp[i,p];
{p in 1..npce[i]} slp[i,p]>> V[i];

To eliminate the product of variables (linearization), I came up with
a scheme which applies only to binary variable 'I[i]' and continuous
variable 'Z' (with lower and upper bound ) following AIMMS ILP tricks,
but in my case since 'Z' variable has only lower bound, thus the
scheme is not appropriate for linearization and thus the results
aren't the ones expected. Below are the definitions of the variables
from my model,

var I {i in NODES} binary; (indicator variable (decision) for
selection of node 'i')
var m {i in NODES, j in NODES} >= 0; (m[i,j]; mapping variable)
var V {i in NODES} = sum {j in NODES} m[i,j] * T [j]; (comment:
T[j]: traffic parameter)
subject to mapp {i in NODES, j in NODES}:
m[i,j] <= I[i]; (constraint: map dst node to source
node only with I[i] set to 1)

Could you provide hints regarding the structure of linearization that
I should follow in my case ?

Thanks!
-- Syed

[Robert Fourer]

You could formulate it as

var Z {i in NODES} >= 0;

minimize max_cost: sum {i in NODES} Z[i];

subject to Zlow {i in NODES}:
Z[i] >= c_f * I[i];

subject to Zref {i in NODES}:
Z[i] >= <<{p in 1..npce[i]-1} bkp[i,p];

{p in 1..npce[i]} slp[i,p]>> V[i] - UB[i] *
(1-I[i]);

where UB[i] is an upper bound on the value that <<{p in 1..npce[i]-1}
bkp[i,p]; {p in 1..npce[i]} slp[i,p]>> V[i] can take in feasible solutions
to your problem. For best chances of getting a good solution within a
reasonable amount of time, choose the smallest upper bound you can think of.

Of course there are no upper bounds in the model as it is currently written,
but in any real problem there is some limit to how large an objective term
can be.

[Syed Hasan]

Hi Bob,
Thank you very much for your answer !
Now, I am getting the expected result from the model : non-linearity
has been solved.
I have set UB[i] to be 100 and 50 (for all nodes) in two different
runs, and the result obtained gave the same objective value (total
cost) in both cases.

Kind Regards,
Syed

On Jun 2, 10:51 pm, "Robert Fourer" <4...@ampl.com> wrote:
> You could formulate it as
>
>    var Z {i in NODES} >= 0;
>    minimize max_cost: sum {i in NODES} Z[i];
>
>    subject to Zlow {i in NODES}:
>                 Z[i] >= c_f * I[i];
>    subject to Zref {i in NODES}:
>                 Z[i] >= <<{p in 1..npce[i]-1} bkp[i,p];
>                          {p in 1..npce[i]} slp[i,p]>> V[i] - UB[i] *
> (1-I[i]);
>
> where UB[i] is an upper bound on the value that <<{p in 1..npce[i]-1}
> bkp[i,p]; {p in 1..npce[i]} slp[i,p]>> V[i] can take in feasible solutions
> to your problem.  For best chances of getting a good solution within a
> reasonable amount of time, choose the smallest upper bound you can think of.
>
> Of course there are no upper bounds in the model as it is currently written,
> but in any real problem there is some limit to how large an objective term
> can be.
>
> Bob Fourer
> 4...@ampl.com
>
> -----Original Message-----
> From: am...@googlegroups.com [mailto:am...@googlegroups.com] On Behalf Of Syed
>
> Hasan
> Sent: Friday, June 01, 2012 4:10 PM
> To: AMPL Modeling Language

> Cc: hasan.syedan...@gmail.com

[Syed Hasan]

Dear Bob, Paul, and all;

Earlier, I have solved my min-max MILP model for 50-80 nodes within feasible time (based on discussion in my earlier thread), but for problem with more than 100 nodes i.e., for problem instance with 200 nodes, the root relaxation phase is taking very long time and getting 'out of memory' issues. I am trying lagrangian relaxation technique to reduce the problem size which will help me to get a feasible solution from the branch and bound algorithm. 

To write the LR algorithm for my model, I have followed both OPL (IBM) and AMPL implementation for the location-transportation problem at links below:

http://www.ensta-paristech.fr/~diam/ro/online/ilog/examples/opl/opl/LagrangianRelaxation/

http://www.ampl.com/NEW/LOOP2/index.html  (trnloc2a.run trnloc2a.mod trnloc2.dat)

and an important document: An Applications Oriented Guide to Lagrangian Relaxation - Marshall Fisher

After following the examples of LR method with the standard subgradient algorithm, I have tried to implement the same for my model, but couldn't formulate properly for the min-max problem, thus getting incorrect results (objective function values are not correct). I know that the lagrangian dual of a minimization problem is a 'max' (lagrangian dual function with multipliers) model, but I do not know how to write an effective and equivalent LR scheme for my model to find/deduce a feasible solution (UpperBound) by employing heuristics. Could you provide me some guidelines/directions on how to apply LR technique for my problem ?.

I am using Gurobi 5.0 solver with AMPL 64-bit.

thanks much!

Syed

#loc_problem.mod

set NODES ordered;  # Vertices

set EDGES within {i in NODES, j in NODES: i != j};

param T {j in NODES};

param D {i in NODES,j in NODES}; # distance-matrix

param mult {NODES} >= 0; # Lagrange multipliers for the constraint (2)

param c_f = 100;

param c_b = 50;

param npce {NODES} integer > 1; # number of pieces for PWLA

param bkp {i in NODES, p in 1..npce[i]-1};  # breakpts and slopes for PWLA

param slp {i in NODES, p in 1..npce[i]};

param UBB {i in NODES};

param d = 1;   

var I {i in NODES} binary;

var m {i in NODES, j in NODES} >= 0;

var V {i in NODES} =  sum {j in NODES} m[i,j] * T [j];

var Z {i in NODES} >= 0;

minimize Lagrangian: sum {i in NODES} Z[i] + sum {j in NODES} mult[j] * (sum {i in NODES} m[i,j] - 1) ;

minimize max_cost: sum {i in NODES} Z[i];

subject to limit {i in NODES}: sum {j in NODES} T[j] * D[i,j] * m[i,j] <= d * V[i];  # performance constraint (1)

subject to start {i in NODES}: sum {j in NODES} m[i,j] = 1;  

# constraint 2: It ensures that a location node should be mapped to each destination node,

# Relaxed for LR relaxation.

                                                                                                              

subject to mapp {i in NODES, j in NODES}:

             m[i,j] <= I[i];

# ensures that destination nodes are mapped only to active location nodes.

# implementing the max (c_f, c_b * V(i)^0.75) cost model:

subject to Zlow {i in NODES}:

                Z[i] >= c_f * I[i];

subject to Zref {i in NODES}:

                Z[i] >= c_b * <<{p in 1..npce[i]-1} bkp[i,p];

                         {p in 1..npce[i]} slp[i,p]>> V[i] - (UBB[i] * (1-I[i]));

Run file:

# Lagrangian relaxation with Subgradient algorithm for location problem

printf "\nLP RELAXATION\n\n";

model loc_problem.mod;

data loc_problem.dat;

option omit_zero_rows 1;

option display_eps .000001;

option solution_round 8;

option solver gurobi_ampl;

option gurobi_options 'lpmethod 3 outlev 1 logfreq 5 threads 5 timing 1';

option log_file '200-nodes-run';

option solver_msg 0; 

option relax_integrality 1;

objective max_cost;

expand;

solve;

param LB; param UB;

let LB := max_cost.val;

let UB := sum {i in NODES} V[i];

option relax_integrality 0;

problem LowerBound: I, m, V, mapp, Lagrangian, limit;

problem UpperBound: m, V, mapp, limit, start, max_cost;

let {j in NODES} mult[j] := 0;

param slack {NODES};

param scale default 1;

param norm;

param step;

param same default 0;

param same_limit := 3;

param iter_limit := 10;

param LBlog {0..iter_limit};  let LBlog[0] := LB;

param UBlog {0..iter_limit};  let UBlog[0] := UB;

param scalelog {1..iter_limit};

param steplog {1..iter_limit};

for {k in 1..iter_limit} { printf "\nITERATION %d\n\n", k;

   solve LowerBound;

   printf "\n solved UB! \n \n";

   printf "\n Lagrangian %f LB %f \n \n", Lagrangian, LB;

   let {j in NODES} slack[j] := sum {i in NODES} m[i,j] - 1;

  if Lagrangian > LB + 0.000001 then {

      let LB := Lagrangian;

      let same := 0; }

   else let same := same + 1;

   if same = same_limit then {

      let scale := scale / 2;

      let same := 0;

      };

   let norm := sum {j in NODES} slack[j]^2;

   printf "\n norm value %f", norm;

   let step := scale * (UB - Lagrangian) / norm;

   let {j in NODES} mult[j] := mult[j] less  step * slack[j];

   if sum {i in NODES} I[i]

         >= 1 - 1e-8 then {

      solve UpperBound;

      let UB := min (UB,max_cost);

      printf "\n solved UB! \n \n";

      printf "\n max_cost %f UB %f \n \n", max_cost, UB;

      }

   let LBlog[k] := LB;

   let UBlog[k] := UB;

   let scalelog[k] := scale;

   let steplog[k] := step;

   }

printf "\n\n";

display LBlog, UBlog, scalelog, steplog;

display I,V,Z;

display max_cost;

display _total_solve_time;