Decision Optimization

 View Only
Expand all | Collapse all

Different results using Cplex in Matlab and Cplex in its own environment

  • 1.  Different results using Cplex in Matlab and Cplex in its own environment

    Posted Fri December 10, 2021 07:08 AM
    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
    ------------------------------

    Attachment(s)



  • 2.  RE: Different results using Cplex in Matlab and Cplex in its own environment

    Posted Fri December 10, 2021 12:02 PM
    It is unclear what you mean by "the power flow is supposed to be negative", but if you are referring to the values of some of the variables (PI, QI?) in the final solution, note that your OPL model declares those variables to be nonnegative (the "+" in "dvar float+").

    ------------------------------
    Paul Rubin
    Professor Emeritus
    Michigan State University
    ------------------------------