Comment appeler le même indice dans chacune des 2 dimensions dans un GAMS variable?

Je suis en cours d'exécution d'un modèle dans GAMS avec un certain nombre de variables avec 2 dimensions, appelant les indices i et j respectivement. J'ai quelques contraintes liées au juste la "diagonale" entrées, donc j'ai codé quelque chose comme "d(i,i)", mais GAMS donne erreur 171 quand je fais cela, par exemple, la première erreur vient sur l'équation e3. Comment puis-je exprimer ce dans GAMBETTES? Code suit:

    Set
   i 'Origin' / 1*20 /
   j 'Destination' / 1*20 /;
   
Table A(i,j)

     1    2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18   19   20
1    0    4    2    574  481  440  408  633  573  1066 486  1926 1537 183  334  374  107  509  378  499 
2    4    0    2    573  480  438  412  632  572  1064 484  1924 1535 182  332  372  111  508  377  498 
3    3    2    0    572  479  437  411  631  571  1063 484  1923 1534 181  332  371  109  507  376  497 
4    574  572  572  0    93   249  135  1188 1128 1620 1040 2480 2091 737  888  928  683  498  367  488 
5    481  480  479  93   0    156  228  1095 1035 1527 948  2387 1998 645  796  836  591  406  275  396 
6    440  438  437  249  156  0    384  1053 993  1485 906  2345 1956 603  754  794  549  364  233  354 
7    409  412  411  135  228  384  0    1044 984  1476 897  2336 1947 594  745  785  302  633  502  623 
8    633  632  631  1188 1095 1053 1044 0    532  1024 623  1246 792  450  300  403  742  560  992  668 
9    573  572  571  1128 1035 993  984  532  0    494  326  879  424  392  385  207  682  1063 932  1053
10   1066 1064 1063 1620 1527 1485 1476 1024 494  0    820  1979 823  885  878  699  1175 1555 1424 1545
11   486  484  484  1040 948  906  897  623  326  820  0    724  270  359  477  223  595  975  844  966 
12   1926 1924 1923 2480 2387 2345 2336 1246 878  1979 724  0    454  1943 1099 846  2035 2415 2284 2405
13   1537 1535 1534 2091 1998 1956 1947 792  424  823  270  454  0    1554 645  392  1646 2026 1895 2016
14   183  182  181  738  645  603  594  450  392  884  359  1943 1554 0    151  192  292  672  542  663 
15   334  332  332  888  796  754  745  300  385  878  477  1099 645  151  0    257  443  261  692  369 
16   373  372  371  928  835  793  784  403  207  699  223  846  392  192  256  0    482  863  732  853 
17   107  111  109  683  591  549  302  742  683  1175 595  2035 1646 292  443  483  0    618  487  608 
18   509  507  507  498  406  364  633  560  1063 1555 975  2415 2026 672  261  863  618  0    134  180 
19   378  377  376  367  275  233  502  992  932  1424 844  2284 1895 542  692  732  487  134  0    122 
20   500  498  497  489  396  354  624  668  1053 1545 966  2405 2016 663  369  854  609  180  122  0   
;

    
Positive variable c(i,j) 'cost direct';
Positive variable cstar(i,j) 'cost routing';
Positive variable z 'objective';
Binary variable d(i,j) 'decision direct';
Binary variable dstar(i,j) 'decision routing';
Binary variable p(i,j) 'd(i,i)*dstar(i,j)';
Binary variable q(i,j) 'd(i,j)*p(i,j)';
Binary variable x(i,j) 'd(j,j)*dstar(i,j)';
Binary variable y(i,j) 'd(i,j)*x(i,j)';
Binary variable f(i,j) 'd(i,i)*d(j,j)';
Binary variable g(i,j) 'f(i,j)*dstar(i,j)';
Integer variable h(i)  'd(i,i)*u(i)';
Integer variable k(i,j) 'd(i,i)*u(j)';
Integer variable m(i,j) 'f(i,j)*u(i)';
Integer variable u(i) 'ordering of tours';

Equation
   e1
   e2(i,j)
   e3
   e4(j)
   e5(i,j)
   e6(i,j)
   e7(i)
   e8(j)
   e9(i)
   e10(i,j)
   e11(i,j)
   e12(i,j)
   e13(i,j)
   e14(i,j)
   e15(i,j)
   e16(i,j)
   e17(i,j)
   e18(i,j)
   e19(i,j)
   e20(i,j)
   e21(i,j)
   e22(i,j)
   e23(i,j)
   e24(i,j)
   e25(i,j)
   e26(i,j)
   e27(i,j)
   e28(i)
   e29(i)
   e30(i)
   e31(i,j)
   e32(i,j)
   e33(i,j)
   e34(i,j)
   e35(i,j)
   e36(i,j)
   e37(i,j)
   e38(i,j)
   e39(i)
   e40(i)
   e41(i)
   e42(i)
   e43(i,j)
   e44(i,j)
   e45(i,j)
   e46(i,j)
   e47(i,j)
   e48(i,j)
   e49(i,j)
   e50(i,j)
   e51(i,j);
   

   e1       .. z =e= sum((i,j),c(i,j))+sum((i,j),cstar(i,j));
   e2(i,j)  .. c(i,j) =e= A(i,j)*d(i,j);
   e3       .. sum((i),d(i,i)) =e= 2;
   e4(j)    .. sum((i),d(i,j)) =e= 1;
   e5(i,j)  .. d(i,j)-d(i,i) =l= 0;
   e6(i,j)  .. cstar(i,j) =e= A(i,j)*dstar(i,j);
   e7(i)    .. dstar(i,i) =e= 0;
   e8(j)    .. sum((i),dstar(i,j)) =e= 1;
   e9(i)    .. sum((j),dstar(i,j)) =e= 1;
   e10(i,j) .. p(i,j) =l= d(i,i);
   e11(i,j) .. p(i,j) =l= dstar(i,j);
   e12(i,j) .. p(i,j) =g= d(i,i)+dstar(i,j)-1;
   e13(i,j) .. p(i,j) =g= 0;
   e14(i,j) .. q(i,j) =l= d(i,j);
   e15(i,j) .. q(i,j) =l= p(i,j);
   e16(i,j) .. q(i,j) =g= d(i,j)+p(i,j)-1;
   e17(i,j) .. q(i,j) =g= 0;
   e18(i,j) .. p(i,j)-q(i,j) =e= 0;
   e19(i,j) .. x(i,j) =l= d(j,j);
   e20(i,j) .. x(i,j) =l= dstar(i,j);
   e21(i,j) .. x(i,j) =g= d(j,j)+dstar(i,j)-1;
   e22(i,j) .. x(i,j) =g= 0;
   e23(i,j) .. y(i,j) =l= d(i,j);
   e24(i,j) .. y(i,j) =l= x(i,j);
   e25(i,j) .. y(i,j) =g= d(i,j)+x(i,j)-1;
   e26(i,j) .. y(i,j) =g= 0;
   e27(i,j) .. x(i,j)-y(i,j) =e= 0;
   e28(i)   .. u(i) =l= 20;
   e29(i)   .. 1000*(1-d(i,i)) =g= u(i)-1;
   e30(i)   .. 2-d(i,i) =l= u(i);
   e31(i,j) .. f(i,j) =l= d(i,i);
   e32(i,j) .. f(i,j) =l= dstar(j,j);
   e33(i,j) .. f(i,j) =g= d(i,i)+dstar(j,j)-1;
   e34(i,j) .. f(i,j) =g= 0;
   e35(i,j) .. g(i,j) =l= f(i,j);
   e36(i,j) .. g(i,j) =l= dstar(i,j);
   e37(i,j) .. g(i,j) =g= f(i,j)+dtar(i,j)-1;
   e38(i,j) .. g(i,j) =g= 0;
   e39(i)   .. 0 =l= h(i);
   e40(i)   .. h(i) =l= 1000*d(i,i);
   e41(i)   .. u(i)-1000*(1-d(i,i)) =l= h(i);
   e42(i)   .. h(i) =l= u(i);
   e43(i,j) .. 0 =l= k(i,j);
   e44(i,j) .. k(i,j) =l= 1000*d(i,i);
   e45(i,j) .. u(j)-1000*(1-d(i,i)) =l= k(i,j);
   e46(i,j) .. k(i,j) =l= u(j);
   e47(i,j) .. 0 =l= m(i,j);
   e48(i,j) .. m(i,j) =l= 1000*f(i,j);
   e49(i,j) .. u(i)-1000*(1-f(i,j)) =l= m(i,j);
   e50(i,j) .. m(i,j) =l= u(i);
   e51(i,j) .. u(i)-u(j)+1-h(i)+k(i,j)-d(i,i)-k(j,i)+h(j)-d(j,j)-m(i,j)+m(j,i)-f(i,j) =l= 19*(1-d(i,i)-d(j,j)-f(i,j)-dstar(i,j)+p(i,j)+x(i,j)+g(i,j));


Model transport / all /;
solve transport using mip minimizing z;
display d.l, dstar.l;