#Constraints
subject to const1{i in V}: sum{j in Pf[i]} f[i,j]-sum{j in Pb[i]} f[j,i]=b1[i]*sum{(entr,j) in E} f[entr,j];
subject to const2{(i,j) in E}: f[i,j]<=q[i,j];
subject to const3{(i,j) in E}: qmin[i,j] <= q[i,j] <= qmax[i,j];
subject to const4{(i,j) in E}: 1 <= 1/x[i,j]*(qmax[i,j]-q[i,j]) <= qmax[i,j]-qmin[i,j];
subject to const5{(i,j) in E}: 1 <=1/y[i,j]* (q[i,j]-qmin[i,j]) <= qmax[i,j]-qmin[i,j];
subject to const6{u in Ag, (i,j) in Eu[u]}: deltaF[i,j,u]= - c[i,j] +(c[i,j]-M)*(1-x[i,j]);
subject to const7{u in Ag, (i,j) in E diff Eu[u]}: deltaF[i,j,u]= -M ;#(i,j) in E but not in Eu[u]
subject to const8{u in Ag, (i,j) in Eu[u]}: deltaB[i,j,u]= c[i,j] -(c[i,j]+M)*(1-y[i,j]) ;
subject to const9{u in Ag, (i,j) in E diff Eu[u]}:deltaB[i,j,u]= -M*(1-y[i,j]);
subject to const10{(i,j) in E,u in Ag}: t[j,u]-t[i,u]>=deltaF[i,j,u];
subject to const11{(i,j) in E,u in Ag}: t[i,u]-t[j,u]>=deltaB[i,j,u];
subject to const12{i in V, u in Ag}: sum{j in Pf[i]} phi[i,j,u]-sum{j in Pb[i]} phi[j,i,u]=b2[i];
subject to const13{u in Ag}: t[exit,u] <= sum{(i,j) in E} psi[i,j,u];
subject to const14{u in Ag}: t[exit,u] <= wu[u] * pi;
subject to const15{(i,j) in E, u in Ag}: psi[i,j,u] >= -M;
subject to const16{(i,j) in E, u in Ag}: psi[i,j,u] <= phi[i,j,u]*c[i,j];
subject to const17{(i,j) in E, u in Ag}: psi[i,j,u] <= deltaF[i,j,u];
subject to const18{(i,j) in E, u in Ag}: psi[i,j,u] >= M*(phi[i,j,u]-1)+deltaF[i,j,u];