Dear Network,

I'm trying to solve a mixed integer linear optimization program for network flow. If I launch the system in CPLEX (using the CPLEX STUDIO 12.10.0 environment) it gives different results than those obtained when using the same problem with exactly the same constraints but using the cplex solver in Matlab.

I've found out that the results are the same when power flow is positive. However, surprisingly, when the power flow is supposed to be negative (source is placed downstream the network) CPLEX does not accept negative values and asign zero to this values. However, when launching in Matlab environment using CPLEX as a solver, it gets the correct results (I attach the results of both cases; CPLEX itself and Matlab using CPLEX solver).

Could it be something related to tolerances? or another internal parameter of the CPLEX solver? This is driving me crazy as getting different results with same constraints and data is something unusual...

Thank you in advance for your help.

Alex

This is the code in CPLEX

- - - - - - -

using CPLEX;

// parameters

int n=5; // number of buses

int m=4; // number of lines

int y=1000; // large number

int z=350; // cost of LS

int t=1; // Rated voltage

float Vmax=1.05;// Max Voltage

float Vmin=0.95;// Min Voltage

range buses = 1..n;// length set of buses

range lines = 1..m;// length set of lines

float PLoad [buses] = [0,0.009,0.009,0.009,0.009];//Max Active Power of dispatchable DG

float QLoad [buses] = [0,0.006,0.006,0.006,0.006];//Max Reactive Power of dispatchable DG

float PDGmax[buses] = [0,0,0,0,0.036];//Active Power OF EACH BUS

float QDGmax[buses] = [0,0,0,0,0.024];//Reactive Power OF EACH BUS

float alpha [lines] = [0,1,1,1];//Alpha coefficient

float r [lines] = [0.093850842,0.025549741,0.044230064,0.044230064];//Lines resistance

float x [lines] = [0.084566834,0.029848586,0.058480517,0.058480517];//Lines reactance

// variables

dvar float+ Pl[lines];//Active power flow across lines

dvar float+ Ql[lines];//Reactive power flow across lines

dvar float+ PDG[buses];//Active power injected by DG

dvar float+ QDG[buses];//Reactive power injected by DG

dvar float+ Pls[buses];//Active powerload shedding flow across line

dvar float+ Qls[buses];//Reactive power load shedding flow across line

dvar float+ V[buses];//Reactive power load shedding flow across line

dvar float+ DELTAV[lines];//Reactive power load shedding flow across line

// model

minimize sum(i in buses) (Pls[i]+Qls[i])*z;

// constraints

subject to

{

forall (i in buses)

0<=Pls[i]<=PLoad[i]; //Limit for active power load shedding

forall (i in buses)

0<=Qls[i]<=QLoad[i]; //Limit for reactive power load shedding

forall (i in buses)

0<=PDG[i]<=PDGmax[i];//Limit for dispatachable DG Active power

forall (i in buses)

0<=QDG[i]<=QDGmax[i]; //Limit for non-dispatachable DG Reactive power

Pl[1] - Pl[2]==(PLoad[2]-PDG[2]-Pls[2]);// ACTIVE power balance at each bus 2

Ql[1] - Ql[2]==(QLoad[2]-QDG[2]-Qls[2]);// ACTIVE power balance at each bus 2

Pl[2] - Pl[3]==(PLoad[3]-PDG[3]-Pls[3]);// ACTIVE power balance at each bus 3

Ql[2] - Ql[3]==(QLoad[3]-QDG[3]-Qls[3]);// ACTIVE power balance at each bus 3

Pl[3] - Pl[4]==(PLoad[4]-PDG[4]-Pls[4]);// ACTIVE power balance at each bus 4

Ql[3] - Ql[4]==(QLoad[4]-QDG[4]-Qls[4]);// ACTIVE power balance at each bus 4

Pl[4]==(PLoad[5]-PDG[5]-Pls[5]);// ACTIVE power balance at each bus 5

Ql[4]==(QLoad[5]-QDG[5]-Qls[5]);// REACTIVE power balance at each bus 5

forall (i in lines)

DELTAV[i]==(Pl[i]*r[i]+Ql[i]*x[i])/t;//Voltage drop across DN lines

forall (i in lines)

(1-alpha[i])*-y<= (V[i]-V[i+1])- DELTAV[i]<=(1-alpha[i])*y;//Voltage drop constraint

forall (j in lines)

-1*alpha[j]<=Pl[j]<=1*alpha[j]; //active power limits

forall (j in lines)

-1*alpha[j]<= Ql[j]<=1*alpha[j];//reactive power limits

forall (i in buses)

Vmin<=V[i]<=Vmax;//Voltage bus limits

V[5]==1;//Imposing voltage of the islanding root DG

}

- - - - - - -

This is the code in MATLAB.

clc;

clear;

%% 1.obtain data

mpc=IEEE4_bus;

Pload=mpc.pload';

Qload=mpc.qload';

PDGmax=mpc.PDGmax';

QDGmax=mpc.QDGmax';

Alpha=mpc.alpha';

r=mpc.r';

x=mpc.x';

%% 2.set parameters

nb=5; %number of buses

nl=4; %number of lines

U=1; %rated voltage

m=350; %cost of load shedding

M=1000; %large number

Vmax=1.05; %maximum voltage

Vmin=0.95; %minimum voltage

%% save results

V_total=zeros(nb,1);

Vl_total=zeros(nl,1);

Pl_total=zeros(nl,1); %active power at each line

Ql_total=zeros(nl,1);

Pls_total=zeros(nb,1);

Qls_total=zeros(nb,1);

PDG_total=zeros(nb,1);

QDG_total=zeros(nb,1);

Cost_total=zeros(1,1);

%% variables

V=sdpvar(nb,1); %voltage at each bus

Vl=sdpvar(nl,1); %voltage drop at each line

Pl=sdpvar(nl,1); %active power at each line

Ql=sdpvar(nl,1); %reactive power at each line

Pls=sdpvar(nb,1); %active power load shedding

Qls=sdpvar(nb,1); %reactive power load shedding

PDG=sdpvar(nb,1); %active power injected by DG

QDG=sdpvar(nb,1); %reactive power injected by DG

%% set constraints

C=[];

%% active power load shedding constraints

for i=1:nb

C=[C,0<=Pls(i),Pls(i)<=Pload(i)];

end

%% reactive power load shedding constraints

for i=1:nb

C=[C,0<=Qls(i),Qls(i)<=Qload(i)];

end

%% DG active power constraints

for i=1:nb

C=[C,0<=PDG(i),PDG(i)<=PDGmax(i)];

end

%% DG reactive power constraints

for i=1:nb

C=[C,0<=QDG(i),QDG(i)<=QDGmax(i)];

end

%% active and reactive power constraints at each bus

C=[C,Pl(1)*Alpha(1)+PDG(2)==Pload(2)-Pls(2)+Pl(2)*Alpha(2)];

C=[C,Ql(1)*Alpha(1)+QDG(2)==Qload(2)-Qls(2)+Ql(2)*Alpha(2)];

C=[C,Pl(2)*Alpha(2)+PDG(3)==Pload(3)-Pls(3)+Pl(3)*Alpha(3)];

C=[C,Ql(2)*Alpha(2)+QDG(3)==Qload(3)-Qls(3)+Ql(3)*Alpha(3)];

C=[C,Pl(3)*Alpha(3)+PDG(4)==Pload(4)-Pls(4)+Pl(4)*Alpha(4)];

C=[C,Ql(3)*Alpha(3)+QDG(4)==Qload(4)-Qls(4)+Ql(4)*Alpha(4)];

C=[C,Pl(4)*Alpha(4)+PDG(5)==Pload(5)-Pls(5)];

C=[C,Ql(4)*Alpha(4)+QDG(5)==Qload(5)-Qls(5)];

%% voltage constraints

for i=1:nb

C=[C,Vmin<=V(i),V(i)<=Vmax];

end

V(5)=1;

%% voltage loss across each line

for l=1:nl

C=[C,Vl(l)==(Pl(l)*r(l)+Ql(l)*x(l))/U];

end

%% voltage loss constraints

for l=1:nl

C=[C,-M*(1-Alpha(l))<=V(l)-V(l+1)-Vl(l),V(l)-V(l+1)-Vl(l)<=M*(1-Alpha(l))];

end

%% active power constraints

for l=1:nl

C=[C,-1*Alpha(l)<=Pl(l),Pl(l)<=Alpha(l)];

end

%% active power constraints

for l=1:nl

C=[C,-1*Alpha(l)<=Ql(l),Ql(l)<=Alpha(l)];

end

%% obkective function

objective=sum(Pls.*m+Qls.*m);

%% solve

ops=sdpsettings('solver','cplex','verbose',0);

sol=optimize(C,objective,ops);

%% save data

Pls_total(:,1)=double(Pls);

Qls_total(:,1)=double(Qls);

Vl_total(:,1)=double(Vl);

Pl_total(:,1)=double(Pl);

Ql_total(:,1)=double(Ql);

PDG_total(:,1)=double(PDG);

QDG_total(:,1)=double(QDG);

V_total=double(V);

Cost_total=double(objective);

%% export data

% filename1='active power of load shedding.xlsx';

% writematrix(Pls_total,filename1);

% filename2='reactive power of load shedding.xlsx';

% writematrix(Qls_total,filename2);

% filename3='cost.xlsx';

% writematrix(Cost_total,filename3);

% filename4='voltage.xlsx';

% writematrix(V_total,filename4);

% filename5='active power flow.xlsx';

% writematrix(Pl_total,filename5);

% filename6='reactive power flow.xlsx';

% writematrix(Ql_total,filename6);

%% analysis error

if sol.problem==0

disp('successful solved');

else

disp('error');

yalmiperror(sol.problem)

end

------------------------------

Alejandro Serrano Fontova

------------------------------