Solving Steiner Tree with GLPK

I'm new at using the GPLK and I'm trying to solve the Steiner Tree problem through it. The mathematical formulation I'm using is this one:

This is the piece of code I'm testing:

# Number of terminal vertexes
param p, integer;

# Number of steiner vertexes
param s, integer;

# Terminal vertexes set
set P := 1..p;

# Steiner vertexes set
set S := (p+1)..(p+s);

# Vertexes set
set V := P union S;

# Edges set
set E, within V cross V;

# Edges cost
param c{(i,j) in E} default 0;

# Variable y
var y{(i,j) in E}, binary;

# Variable x
var x{(i,j) in E} >= 0;

# x constraint
s.t. Xconstraint{(i,j) in E}: x[i,j] <= (p) * y[i,j];

# Steiner vertexes constraints
s.t. Sconstraint{i in S}: sum{(i,j) in E} x[i,j] - sum{(k,i) in E} x[k,i] = 0;

# Terminal vertexes constraints
s.t. Pconstraint{i in P}: sum{(i,j) in E} x[i,j] - sum{(k,i) in E} x[k,i] = -1;

# Minimize vertexes cost
minimize cost: sum{(i,j) in E} c[i,j]*y[i,j];

solve;

printf "The tree cost is:\n",
   sum{(i,j) in E} c[i,j] * y[i,j];
printf "\n";

data;

param p := 10;
param s := 8;

param : E : c :=
    1 16 1340771
    2 18 2783118
    2 14 1534253
    2 17 71057
    3 12 1439171
    3 15 1921785
    4 11 3793393
    4 13 573690
    5 18 1268333
    6 17 2907508
    7 11 1981430
    8 11 1205285
    8 12 3190105
    9 12 393839
    10 16 461534
    10 18 1072719
    11 4 3793393
    11 7 1981430
    11 8 1205285
    12 3 1439171
    13 4 573690
    12 8 3190105
    12 9 393839
    13 14 776147
    13 17 2239343
    14 2 1534253
    14 13 776147
    14 15 2613116
    15 3 1921785
    15 14 2613116
    15 16 170002
    16 1 1340771
    16 10 461534
    16 15 170002
    17 2 71057
    17 6 2907508
    17 13 2239343
    18 2 2783118
    18 5 1268333
    18 10 1072719
;

end;

I've tested this exemple using this site: http://www3.nd.edu/~jeff/mathprog/mathprog.html and It display this message: "PROBLEM HAS NO FEASIBLE SOLUTION". I believe either this example is bad or the constraints are wrong, can someone help me?

You need to change this constraint:

# Terminal vertexes constraints
s.t. Pconstraint{i in P: i != 1}: sum{(i,j) in E} x[i,j] - sum{(k,i) in E} x[k,i] = -1;

That's because the first node is where flows from in this flow based formulation.