Feasible solution opens arcs which have no flow over them

[b.g.van...@gmail.com]

Hi, 

When I solve the following model using KNITRO, a feasible solution is found within the specified timelimit. 

However, when I check the solutions, I noticed that when there is no flow x[i,j,k] over certain arcs (i,j), the arcs are still opened (indicated by y[i,j]) and fixed costs y[i,j] * f[i,j] are incurred. 

This seems very strange to me, since the objective is to minimize costs. Furthermore, I solved the model with a linear decay function in the objective instead of an exponential decay function but with the exact same constraints using CPLEX, and in this only arcs are opened which are in use. 

Is this a problem due to the solver being used, or am I missing any extra constraints which should be included when using an exponential instead of a linear cost function? Thanks in advance

.mod file with exponential cost function:

# CMND WITH LIFETIME RESTRICTIONS AND AN EXPONENTIAL QUALITY LOSS FUNCTION #

### Sets ###

set N; # set of nodes

set A within N cross N; # set of directed arcs

set S within N; # set of suppliers

set D within N; # set of demand nodes

set T:= {i in N: i not in S and i not in D}; # set of transshipment points

set P; # set of products

set K; # set of commodities representing due dates

### Parameters ###

param xcoord{i in N}; # x-coordinate of node i

param ycoord{i in N}; # y-coordinate of node i

param f{(i,j) in A} >=0; # fixed cost applied on arc (i,j)

param u{(i,j) in A} >=0; # capacity of arc (i,j)

param supply{i in S, k in K} >=0; # supply of commodity k from supplier s

param demand{i in D, p in P} >=0; # demand in node i for product p

param distance{(i,j) in A}:= round((sqrt((xcoord[i] - xcoord[j])^2 + 

  (ycoord[i] - ycoord[j])^2))); # distance over arc (i,j)

param t{(i,j) in A}:= round(distance[i,j]/1.2); # travel time over arc (i,j) in hours

param tc:= 0.0179; # transportation cost per unit per hour

param rcf; # relative cost factor

param c{(i,j) in A, k in K}:= tc * t[i,j] * rcf; # variable transportation cost per unit over arc (i,j)

param E{k in K} symbolic in P; # the product represented by commodity k 

param b{(i,j) in A, p in P}:= min(demand[i,p],u[i,j]); # an upper bound of the flow of product p on arc (i,j)

param slmin{k in K} >=0; # point in time where commodity k starts deteriorating in days

param slmax{k in K} >=0; # maximum shelf life of commodity k in days

param M; # "Big M": sufficiently large number

param price{k in K} >=0; # price of commodity k

param rate{k in K}; # decay rate per hour in dollars

### Decision variables ###

var x{(i,j) in A, k in K} >=0 integer; # flow distribution decision variable indicating the amount of flow of commodity k on arc (i,j)

var y{(i,j) in A} binary; # design variable; equals 1 if (i,j) is selected in the final design and 0 if otherwise

var z{(i,j) in A, k in K} binary; # 1 if commodity k travels on arc (i,j), 0 if otherwise

var a{i in N, k in K} >=0; # arrival time of commodity k in node i in days

var w{i in N, k in K} >=0; # arrival time of commodity k in destination node i in days

### Objective function ###

minimize TotalCost:

sum{(i,j) in A, k in K} c[i,j,k] * x[i,j,k] + # transportation costs

sum{(i,j) in A} f[i,j] * y[i,j] + # fixed costs

sum{k in K, i in D, p in P: p = E[k]} demand[i,p] * price[k] * (1 - exp(-rate[k] * w[i,k])) + # decay costs

sum{i in N, k in K} 0.0001 * w[i,k] +

sum{i in N, k in K} 0.0001 * a[i,k];

### Constraints ###

s.t. DemandConstraint{p in P, i in D}:

sum{k in K: E[k] = p} (sum{(j,i) in A} x[j,i,k] - sum{(i,j) in A} x[i,j,k]) = demand[i,p];

s.t. SupplyConstraint{k in K, i in S}:

sum{(i,j) in A} x[i,j,k] - sum{(j,i) in A} x[j,i,k] <= supply[i,k];

s.t. FlowConservationConstraint{k in K, i in T}:

sum{(i,j) in A} x[i,j,k] - sum{(j,i) in A} x[j,i,k] = 0;

s.t. ArcCapacity{(i,j) in A}:

sum{k in K} x[i,j,k] <= u[i,j] * y[i,j];

s.t. ArrivalTimeForEachPath{k in K, (i,j) in A}:

a[i,k] + t[i,j] <= a[j,k] + M * (1 - z[i,j,k]);

s.t. LifetimeRestriction{i in N, k in K}:

a[i,k] <= slmax[k];

s.t. DecayPenalty{i in D, k in K}:

a[i,k] <= slmin[k] + w[i,k];

s.t. MaxFlow1{k in K, (i,j) in A}:

x[i,j,k] <= M * z[i,j,k];

s.t. MaxFlow2{k in K, (i,j) in A}:

z[i,j,k] <= y[i,j];

s.t. Starttime{i in S, k in K}:

a[i,k] = 0;

.mod file with linear cost function:

# CMND WITH LIFETIME RESTRICTIONS AND A LINEAR LOSS FUNCTION #

### Sets ###

set N; # set of nodes

set A within N cross N; # set of directed arcs

set S within N; # set of suppliers

set D within N; # set of demand nodes

set T:= {i in N: i not in S and i not in D}; # set of transshipment points

set P; # set of products

set K; # set of commodities representing due dates

### Parameters ###

param xcoord{i in N}; # x-coordinate of node i

param ycoord{i in N}; # y-coordinate of node i

param f{(i,j) in A}:= round(Uniform(1000,3000)) >=0; # fixed cost applied on arc (i,j)

param u{(i,j) in A}:= round(Uniform(10000,50000)) >=0; # capacity of arc (i,j)

param supply{i in S, k in K}:= round(Uniform(0,10000)) >=0; # supply of commodity k from supplier s

param demand{i in D, p in P}:= round(Uniform(2500,7500)) >=0; # demand in node i for product p

param tc:= 0.0179:= Uniform(0.41,0.50); # transportation cost per unit per day

param distance{(i,j) in A}:= round((sqrt((xcoord[i] - xcoord[j])^2 + 

  (ycoord[i] - ycoord[j])^2))); # distance over arc (i,j)

param t{(i,j) in A}:= round(distance[i,j])/1.2; # travel time over arc (i,j)

param rcf; # relative cost factor

param c{(i,j) in A, k in K}:= tc * t[i,j] * rcf; # variable transportation cost per unit over arc (i,j)

param E{k in K}; # the product represented by commodity k 

param slmin{k in K} >=0; # point in time where commodity k starts deteriorating in days

param slmax{k in K} >=0; # maximum shelf life of commodity k in days

param M; # "Big M": sufficiently large number

param price{k in K} >=0; # price of commodity k

param rate{k in K} >=0; # penalty rate per day

### Decision variables ###

var x{(i,j) in A, k in K} >=0 integer; # flow distribution decision variable indicating the amount of flow of commodity k on arc (i,j)

var y{(i,j) in A} binary; # design variable; equals 1 if (i,j) is selected in the final design and 0 if otherwise

var z{(i,j) in A, k in K} binary; # 1 if commodity k travels on arc (i,j), 0 if otherwise

var a{i in N, k in K} >=0; # arrival time of commodity k in node i in days

var w{i in N, k in K} >=0; # arrival time of commodity k in destination node i in days

### Objective function ###

minimize TotalCost:

sum{(i,j) in A, k in K} c[i,j,k] * x[i,j,k] + # transportation costs

sum{(i,j) in A} f[i,j] * y[i,j] + # fixed costs

sum{k in K, i in D, p in P: p = E[k]} (rate[k] * demand[i,p] * price[k] * w[i,k]) + # decay costs

sum{i in N, k in K} 0.0001 * w[i,k];

### Constraints ###

s.t. DemandConstraint{p in P, i in D}:

sum{k in K: E[k] = p} (sum{(j,i) in A} x[j,i,k] - sum{(i,j) in A} x[i,j,k]) = demand[i,p];

s.t. SupplyConstraint{k in K, i in S}:

sum{(i,j) in A} x[i,j,k] - sum{(j,i) in A} x[j,i,k] <= supply[i,k];

s.t. FlowConservationConstraint{k in K, i in T}:

sum{(i,j) in A} x[i,j,k] - sum{(j,i) in A} x[j,i,k] = 0;

s.t. ArcCapacity1{(i,j) in A}:

sum{k in K} x[i,j,k] <= u[i,j] * y[i,j];

s.t. ArrivalTimeForEachPath{k in K, (i,j) in A}:

a[i,k] + t[i,j] <= a[j,k] + M * (1 - z[i,j,k]);

s.t. LifetimeRestriction{i in N, k in K}:

a[i,k] <= slmax[k];

s.t. DecayPenalty{i in D, k in K}:

a[i,k] <= slmin[k] + w[i,k];

s.t. MaxFlow1{k in K, (i,j) in A}:

x[i,j,k] <= M * z[i,j,k];

s.t. MaxFlow2{k in K, (i,j) in A}:

z[i,j,k] <= y[i,j];

s.t. Starttime{i in S, k in K}:

a[i,k] = 0;

[b.g.van...@gmail.com]

Additionally, when I solve the model for a long time, less arcs are being opened but still some unused arcs remain open. Is this a case of Knitro settling for a feasible point due to the timelimit, or am I missing an important constraint? I can't identify which constraint should be added or modified though, since the model with linear decay works just fine.

[BVG]

I also tried to add the following constraint

s.t. MaxFlow3{k in K, (i,j) in A}:

y[i,j] <= x[i,j,k];

In addition to the other flow constraints:

s.t. ArcCapacity1{(i,j) in A}:

sum{k in K} x[i,j,k] <= u[i,j] * y[i,j];

s.t. MaxFlow1{k in K, (i,j) in A}:

x[i,j,k] <= M * z[i,j,k];

s.t. MaxFlow2{k in K, (i,j) in A}:

z[i,j,k] <= y[i,j];

From my understanding, this should ensure that no unused arcs are opened when no flow is directed over that arc. Would this be a valid constraint or am I missing something here?

[Robert Fourer]

Be sure that Knitro reported "optimal solution" since otherwise it is not claiming the solution to be optimal. Then to check whether there is an obviously better solution, after solving, choose some y-variable that is 1 but that you think should be 0, an set it so zero; for example,

let y['A','B',3] := 0;

Then check the slacks for feasibility as explained previously in this group, and use "display TotalCost;" to check whether the objective value has decreased.

If Knitro claimed optimality, and the problem is still feasible and the objective has decreased, then possibly Knitro found a locally optimal solution that is not globally optimal. This can happen when the objective function is not convex. However for more advice you could consider sending a complete example to the developers of Knitro at http://groups.google.com/group/knitro.

Bob Fourer
am...@googlegroups.com

=======

[BVG]

I dropped some of the arcs as you suggested in the following way:

let y['T2','D6'] := 0;

And it gave an improved total cost. However, when I tried to let more binaries equal 0, e.g.:

let y['T2','D6'] := 0;
let y['T1','D3'] := 0;
let y['S1','D4'] := 0;

The feasible solution still opens these arcs:

ampl: display y;

y [*,*] (tr)

:   S1  S2  S3  S4  S5  T1  T2  T3    :=

D1   .   .   .   .   .   1   0   0

D2   .   .   .   .   .   1   0   0

D3   .   .   .   .   .   1   0   0

D4   .   .   .   .   .   0   1   0

D5   .   .   .   .   .   0   0   1

D6   .   .   .   .   .   0   1   1

T1   1   0   1   1   1   .   .   .

T2   0   1   0   1   0   .   .   .

T3   1   0   0   1   0   .   .   .

;

Furthermore you suggested typing the following command:

let y['A','B',3] := 0;

However, my y[i,j] has only 2 subscripts, so I'm not sure what the '3' should do?

Thanks in advance. 

Op zondag 16 juli 2017 15:13:08 UTC+2 schreef Robert Fourer:

[Robert Fourer]

Sorry, I meant "let y['T2','D6'] := 0;" -- just two subscripts. This is just a test to determine whether Knitro claimed optimality but there's actually a better solution. You need to solve once, check that Knitro claimed "optimal solution", then change one y-value from 1 to 0, then check feasibility and the objective value of the changed solution (without solving again).

Bob Fourer
am...@googlegroups.com

[BVG]

I have not yet been able to solve the model to optimality, only to feasible solutions. I have run the problem a couple of times for 10 hours and either the timelimit is reached or the solver runs out of memory. 

Furthermore I'm not sure how to check the feasibility and objective value with  "let y['T2','D6'] := 0;", without solving it again.

[BVG]

In addition to my previous post, would the constraint;

s.t. maxflow3{k in K, (i,j) in A}:

y[i,j] <= x[i,j,k];

be a valid constraint to ensure no unused arcs are opened?

[Robert Fourer]

A Knitro run may end because the time limit has been reached, or it may end because it has determined (local) optimality of the current solution. You can tell which termination condition was reached by the message that you see at the end of the Knitro run:

Knitro 10.2.0: Locally optimal solution.

or

Knitro 10.2.0: Time limit reached. Current point is feasible.

If you get the "Time limit reached" message then you should not be surprised if the solution returned is not optimal. If you get the "Locally optimal solution" message then you can be confident that the solution is optimal among all nearby solutions, but if your objective function is not convex (and in fact, yours appears to be concave) then there is still the possibility that better solutions exist. In that case here are two possibilities:

-- You could add constraints that "x[i,j,k] = 0 for all k" implies "y[i,j] = 0". Since the x-variables are integer and the y-variables are binary, this can be written "y[i,j] <= sum {k in K} x[i,j,k]".

-- You could try BARON or Couenne which are solvers designed to find global solutions. However global optimality is a much harder problem and you would have to run some tests to determine whether you could get optimal solutions from these solvers in any reasonable amount of time.

Bob Fourer
am...@googlegroups.com

-----Original Message-----
From: am...@googlegroups.com [mailto:am...@googlegroups.com] On Behalf Of b.g.van...@gmail.com
Sent: Sunday, July 16, 2017 2:24 AM
To: AMPL Modeling Language

[BVG]

Thanks for your help. 

I have run the model many times for over 8 hours with Knitro, but it has never reached an optimal solution. 

I have tried the BARON solver and although I have not been able to find optimal solutions in a reasonable amount of time, I get much better results than with other solvers, even at relatively low timelimits. 

[BVG]

In addition to my previous post. I have run the model for over 10 hours and this gives the same results as running it for 15 minutes. 

Is there any way to know whether much better solutions exist, or can it be assumed that if there is no change in the results for such a long time, the value has stabilized?

[Robert Fourer]

There are ways of knowing that the best solution has been found, but they are already implemented in the solvers. Thus if the solver continues to run, it's because the current solution cannot yet be proved to be optimal by the methods that the solver uses. It is not unusual that the time to prove optimality is very much longer than the time to find the best solution. (And in the case of LocalSolver, the methods used may not be able to construct any proof of optimality.)

Thus when the solver does not find any better solutions for a long time, you may reasonably stop the solver and use the best solution that has been found. Although optimality has not been proved, there is a good chance that the solution is optimal, or at least very good.

Bob Fourer
am...@googlegroups.com

=======

From: am...@googlegroups.com [mailto:am...@googlegroups.com] On Behalf Of BVG
Sent: Monday, July 24, 2017 10:15 AM
To: AMPL Modeling Language