Drone and truck delivery problem

Hello,
I am writing a mixed integer problem in AMPL, trying to develop a model for a problem of trucks+drones delivery found in literature (the first proposed in the attached paper). Unfortunately, I’m having problems as AMPL gives me errors that I do not completely understand, especially it says that ‘all variables have been eliminated’ in some constraints. You can find attached my .mod and .dat files.
Any suggestions and help would be really appreciated.
Thank you very much in advance,
Lara

FILE.MOD
set nodes;

set D; #depot

set U; #parkings

set W; #customers

set K; #trucks

param coord_x {nodes};

param coord_y {nodes};

param b {nodes}; #demand

param theta;

param Q; #trucks capacity

set E := {i in nodes, j in nodes: i != j};

set A:= {s in W, i in U};

param dist {nodes, nodes};

param dist1 {W, U};

param p {nodes} default 0;

param c {nodes, nodes};

param d {nodes, nodes};

param R;

var z {i in nodes} binary;

var x {i in nodes, j in nodes, k in K : i != j} binary; #: i != j

var y {s in W, i in U} binary;

var q {i in nodes, k in K} >=0;

minimize obj : sum {k in K} sum {(i,j) in E} c[i,j] * x[i,j,k] + sum {s in W, i in U} d[s,i]*y[s,i] + sum {i in U} p[i]*z[i];

subject to constr_b : z[1] = 1;

subject to constr_c {k in K} : sum { (i,j) in E : j != 1} x[1,j,k] <= 1;

subject to constr_d {i in nodes} : sum {k in K} sum {j in nodes: j != i} x[i,j,k] = z[i];

subject to constr_e {s in W} : sum {i in U} y [s,i] = 1 - z[s];

subject to constr_f {i in nodes, k in K} : sum {j in nodes: j != i} x[i,j,k] = sum {j in nodes: j != i} x[j,i,k];

subject to constr_g {s in W, i in U} : d[s,i] * y [s,i] <= R * z[i];

#subject to constr_h {s in W} : y [s,1] = 0;

subject to constr_i {i in nodes, j in W, k in K : i != j} : q[j,k] >= q[i,k] + b[j] + Q * (x[i,j,k]-1);

subject to constr_j {i in nodes, j in U, k in K: i != j} : q[j,k] >= q[i,k] + sum {s in W} b[s] * y [s,j] + Q * ( x[i,j,k] - 1);

subject to constr_k {i in nodes, k in K} : q[i,k] <= Q * z[i];

FILE.DAT
set nodes:= 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26;

set D:= 1;

set U:= 2 3 4 5 6 7;

set W:= 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26;

set K:= 1 2 3;

param coord_x:=

1 30

2 37

3 49

4 52

5 20

6 40

7 21

8 17

9 31

10 52

11 51

12 42

13 31

14 5

15 12

16 36

17 52

18 27

19 17

20 13

21 57

22 62

23 42

24 16

25 8

26 7;

param coord_y:=

1 40

2 52

3 49

4 64

5 26

6 30

7 47

8 63

9 62

10 33

11 21

12 41

13 32

14 25

15 42

16 16

17 41

18 23

19 33

20 13

21 58

22 42

23 57

24 57

25 52

26 38;

param b:=

1 0

2 7

3 30

4 16

5 9

6 21

7 15

8 19

9 23

10 11

11 5

12 19

13 29

14 23

15 21

16 10

17 15

18 3

19 41

20 9

21 28

22 8

23 8

24 16

25 10

26 28;

param theta:= 5;

param Q:= 35;

for {i in nodes, j in nodes}{

let dist[i,j]:= sqrt((coord_x[i] - coord_x[j])^2 + (coord_y[i]-coord_y[j])^2);

let c[i,j]:= dist[i,j]*theta;

}

for {s in W, i in U}{

let dist1[s,i]:= sqrt((coord_x[s] - coord_x[i])^2 + (coord_y[s]-coord_y[i])^2);

let d[s,i]:= (10 - theta) * dist1[s,i];

}

let R:= max {s in W, i in U} dist1[s,i];


By analyzing the constraints of your problem, AMPL’s presolve phase can often determine tighter bounds for the variables. Sometimes presolve determines that a variable’s lower bound and its upper bound are equal; in that case, presolve can fix the variable at the determined bound (and can eliminate the variable from the problem). In your case, it happened that some constraints had all of their variables fixed — and they were fixed at values that did not satisfy the constraint! This proves that no feasible solution to your constraints is possible, and as a result, AMPL did not send your problem to the solver.

As a first step in determining why no feasible solution is possible, you can use AMPL’s expand command to see the constraints that AMPL actually generated. Here are some examples:

ampl: expand constr_b;
subject to constr_b:
	z[1] = 1;

ampl: expand constr_c[1];
subject to constr_c[1]:
	25*x[1,2,1] + 25*x[1,3,1] + 25*x[1,4,1] + 25*x[1,5,1] + 25*x[1,6,1] + 
	25*x[1,7,1] + 25*x[1,8,1] + 25*x[1,9,1] + 25*x[1,10,1] + 25*x[1,11,1]
	 + 25*x[1,12,1] + 25*x[1,13,1] + 25*x[1,14,1] + 25*x[1,15,1] + 
	25*x[1,16,1] + 25*x[1,17,1] + 25*x[1,18,1] + 25*x[1,19,1] + 
	25*x[1,20,1] + 25*x[1,21,1] + 25*x[1,22,1] + 25*x[1,23,1] + 
	25*x[1,24,1] + 25*x[1,25,1] + 25*x[1,26,1] <= 1;

ampl: expand constr_d[1];
subject to constr_d[1]:
	-z[1] + x[1,2,1] + x[1,2,2] + x[1,2,3] + x[1,3,1] + x[1,3,2] + x[1,3,3]
	 + x[1,4,1] + x[1,4,2] + x[1,4,3] + x[1,5,1] + x[1,5,2] + x[1,5,3] + 
	x[1,6,1] + x[1,6,2] + x[1,6,3] + x[1,7,1] + x[1,7,2] + x[1,7,3] + 
	x[1,8,1] + x[1,8,2] + x[1,8,3] + x[1,9,1] + x[1,9,2] + x[1,9,3] + 
	x[1,10,1] + x[1,10,2] + x[1,10,3] + x[1,11,1] + x[1,11,2] + x[1,11,3]
	 + x[1,12,1] + x[1,12,2] + x[1,12,3] + x[1,13,1] + x[1,13,2] + 
	x[1,13,3] + x[1,14,1] + x[1,14,2] + x[1,14,3] + x[1,15,1] + x[1,15,2]
	 + x[1,15,3] + x[1,16,1] + x[1,16,2] + x[1,16,3] + x[1,17,1] + 
	x[1,17,2] + x[1,17,3] + x[1,18,1] + x[1,18,2] + x[1,18,3] + x[1,19,1]
	 + x[1,19,2] + x[1,19,3] + x[1,20,1] + x[1,20,2] + x[1,20,3] + 
	x[1,21,1] + x[1,21,2] + x[1,21,3] + x[1,22,1] + x[1,22,2] + x[1,22,3]
	 + x[1,23,1] + x[1,23,2] + x[1,23,3] + x[1,24,1] + x[1,24,2] + 
	x[1,24,3] + x[1,25,1] + x[1,25,2] + x[1,25,3] + x[1,26,1] + x[1,26,2]
	 + x[1,26,3] = 0;

ampl: expand constr_e[8];
subject to constr_e[8]:
	z[8] + y[8,2] + y[8,3] + y[8,4] + y[8,5] + y[8,6] + y[8,7] = 1;

You can examine these constraints to see whether they are what you expect. If you find that one of these constraints is not right, then by examining it carefully, you can determine what you need to fix in the corresponding subject to statement of the AMPL model.

(You can also use a command like expand constr_d; to see all of the constraints generated from one constraint definition, or expand constr_d >constr_d.txt; to write all the constraints to a file.)

Also some solvers have features to help you find the cause of the infeasibility, but get help with that, you need to say which solver you are using.