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
------------------------------
#DecisionOptimization