Because I am unable to do it in a single program, the values of pevch [j] [k], pevdis [j] [k], pbat [j], and c [j] are actually input from the outputs of another optimization program. This program's decision variables are pbatch[j][k], pbatnew[j], and so on.
int nb=14;
int t=24;
int nv=400;
int l=14;
range nbus=1..nb;
range tavail =1..t;
range nport=1..nv;
//soc data
float soc[nport][tavail]=...;
//forcasted load data
float pl[nbus][tavail]=...;
//target load data
float pt[nbus][tavail]=...;
//generation
float pg[nbus][tavail]=...;
//price data
float c[tavail]=[3.5,3.5,2.5,2.5,2.5,2.5,2.5,2.5,3.5,3.5,5,5,5,5,5,5,5,5,3.5,5,5,5,3.5,3.5];
//decision variable charging power
dvar float+ pbat[tavail];
dvar float+ pbatnew[tavail];
dvar float+ pbatch[tavail][nport];
float G[nbus][nbus]=...;
float B[nbus][nbus]=...;
float deltacos[nbus][nbus][tavail];
float deltasin[nbus][nbus][tavail];
float v[tavail][nbus]=...;
float del[tavail][nbus]=...;
float p[nbus][tavail]=...;
float pevch[tavail][nport]=...;
float pevdis[tavail][nport]=...;
execute fill_deltacos
{for(var i in nbus)for (var l in nbus)for(var j in tavail)deltacos[i][l][j]=Math.cos(del[j][i]-del[j][l]);
}
execute fill_deltasin
{for(var i in nbus)for (var l in nbus)for(var j in tavail)deltacos[i][l][j]=Math.sin(del[j][i]-del[j][l]);
}
minimize 0.7*sum(i in nbus)(sum(j in tavail)sum(k in nport)((pt[i][j]*1000000-pl[i][j]*1000000)+sum(k in nport)(pevdis[j][k]-pevch[j][k]-pbatch[j][k])))- 0.15*(sum (j in tavail)(sum (k in nport)(c[j]*pevdis[j][k])))
+0.15* (sum (j in tavail)sum (k in nport)(c[j]*pevch[j][k]+3*pbatch[j][k]));
subject to
{
//constraint for discharging
forall(i in nbus,j in tavail,k in nport)
{
pbat[j]==sum (i in nbus) p[i][j]*1000000+ sum(i in nbus)(pl[i][j]-pg[i][j])*1000000+sum (k in nport)(-pevch[j][k]+pevdis[j][k]);
ct1:if(c[j]==5&&soc[k][j]<0.7)
{pbatch[j][k]==(0.7-soc[k][j])*100000;
pbatnew[j]==pbat[j]-(sum (k in nport)pbatch[j][k]);}
}
}
