Problem with the OF when replacing =g= with =e=

moatta · October 14, 2019, 9:31pm

Dear Friends,
I have a problem with this code. My objective is to minimize the cost in order to find the optimal power flow in a power system. The cost is equal to the output power from all generators multiplied by the unit cost of each generator. So, it should be =e=. However, when I put =e=, the solution is 9 which is not logical at all because the least MW costs 10 and I have 10 MWs. When I put =g=, the solution becomes 1056 and it is logical. So, what is the problem?

Set
bus / 13 /
slack(bus) / 3 /
Gen / g1g3 /;

Alias (bus,node);

Table GenData(Gen,*) ‘generating units characteristics’
b pmin pmax
g1 10 0 65
g2 11 0 100;

Set GBconect(bus,Gen) ‘connectivity index of each generating unit to each bus’ / 1.g1, 3.g2 /;

Table BusData(bus,*) ‘demands of each bus in MW’
Pd
2 100;

Set conex ‘bus connectivity matrix’ / 1.2, 2.3, 1.3 /;
conex(bus,node)$(conex(node,bus)) = 1;

Table branch(bus,node,*) ‘network technical characteristics’
x Limit
1.2 0.2 50
2.3 0.25 100
1.3 0.4 100 ;

branch(bus,node,‘x’)(branch(bus,node,'x')=0) = branch(node,bus,'x'); branch(bus,node,'Limit')(branch(bus,node,‘Limit’)=0) = branch(node,bus,‘Limit’);
branch(bus,node,‘bij’)$conex(bus,node) = 1/branch(bus,node,‘x’);

Variable OF, Pij(bus,node), Pg(Gen), delta(bus);
Equation const1, const2, const3, const9, const10, const4, const5, const6, const7, const8;

const1(bus,node)$conex(bus,node)…
Pij(bus,node) =e= branch(bus,node,‘bij’)*(delta(bus) - delta(node));

const2(bus)…
sum(Gen$GBconect(bus,Gen), Pg(Gen)) - BusData(bus,‘Pd’) =e= sum(node$conex(node,bus), Pij(bus,node));

const3…
OF =e= sum(Gen,Pg(Gen)*GenData(Gen,‘b’));

const9(bus,node)$conex(bus,node)…
Pij(bus,node)=l=branch(bus,node,‘Limit’);

const10(bus,node)$conex(bus,node)…
Pij(bus,node)=g=-branch(bus,node,‘Limit’);

const4(Gen)…
Pg(Gen)=l=GenData(Gen,‘Pmax’);

const5(Gen)…
Pg(Gen)=g=GenData(Gen,‘Pmin’);

const6(bus)…
delta(bus)=l=pi;

const7(bus)…
delta(bus)=g=-pi;

const8(bus)…
delta(‘3’)=e=0;

Model loadflow / all /;

solve loadflow minimizing OF using lp;
display Pij.l, of.l, pg.l, delta.l;

Renger · October 15, 2019, 6:41am

Hi
If I run your model with the equal sign, I get an infeasible solution as the constraint 7 is infeasible. Taking this constraint out leads to an optimal solution.
Set =E= to =G= doesn’t change the problem of infeasibility.
So I don’t know if you are running the same code.
If you post code, put it between

and

so it is possible to keep your code correctly formatted (especially the tables).

Cheers
Renger

moatta · October 15, 2019, 7:09am

Hello Renger,
Thanks for your reply. I am running the same code and I am getting solutions. Here is the code.

Set
   bus        / 1*3   /
   slack(bus) / 3     /
   Gen        / g1*g3 /;

Alias (bus,node);

Table GenData(Gen,*) 'generating units characteristics'
       b   pmin  pmax
   g1  10  0     65
   g2  11  0     100;
* -----------------------------------------------------

Set GBconect(bus,Gen) 'connectivity index of each generating unit to each bus' / 1.g1, 3.g2 /;

Table BusData(bus,*) 'demands of each bus in MW'
       Pd
   2  100;

Set conex 'bus connectivity matrix' / 1.2, 2.3, 1.3 /;
conex(bus,node)$(conex(node,bus)) = 1;

Table branch(bus,node,*) 'network technical characteristics'
        x     Limit
   1.2  0.2   50
   2.3  0.25  100
   1.3  0.4   100 ;

branch(bus,node,'x')$(branch(bus,node,'x')=0)         =   branch(node,bus,'x');
branch(bus,node,'Limit')$(branch(bus,node,'Limit')=0) =   branch(node,bus,'Limit');
branch(bus,node,'bij')$conex(bus,node)                = 1/branch(bus,node,'x');
*****************************************************

Variable OF, Pij(bus,node), Pg(Gen), delta(bus);
Equation const1, const2, const3, const9, const10, const4, const5, const6, const7, const8;


const1(bus,node)$conex(bus,node)..
   Pij(bus,node) =e= branch(bus,node,'bij')*(delta(bus) - delta(node));

const2(bus)..
   sum(Gen$GBconect(bus,Gen), Pg(Gen)) - BusData(bus,'Pd') =e= sum(node$conex(node,bus), Pij(bus,node));

const3..
   OF =e= sum(Gen,Pg(Gen)*GenData(Gen,'b'));

const9(bus,node)$conex(bus,node)..
         Pij(bus,node)=l=branch(bus,node,'Limit');

const10(bus,node)$conex(bus,node)..
         Pij(bus,node)=g=-branch(bus,node,'Limit');

const4(Gen)..
         Pg(Gen)=l=GenData(Gen,'Pmax');

const5(Gen)..
         Pg(Gen)=g=GenData(Gen,'Pmin');

const6(bus)..
         delta(bus)=l=pi;

const7(bus)..
         delta(bus)=g=-pi;

const8(bus)..
         delta('3')=e=0;

Model loadflow / all /;

solve loadflow minimizing OF using lp;
display Pij.l, of.l, pg.l, delta.l;

Renger · October 16, 2019, 11:24am

Hi
If I run it with CPLEX, I get an infeasible model:

Iteration log . . .
Iteration:     1    Infeasibility =            16.858406
LP status(3): infeasible
Cplex Time: 0.00sec (det. 0.02 ticks)

Model has been proven infeasible.
--- Restarting execution
--- Reading solution for model loadflow
*** Status: Normal completion
--- Job cge2.gms Stop 10/16/19 13:21:47 elapsed 0:00:03.647
 
GAMS process finished at Wed Oct 16 13:21:49 2019
Total running time is 00:00:08.

With constraint const7 as the culprit

---- EQU const7  

     LOWER     LEVEL     UPPER    MARGINAL

1    -3.142    -2.500     +INF       .         
2    -3.142   -12.500     +INF      1.000 INFES
3    -3.142      .        +INF       .         

**** REPORT SUMMARY :        0     NONOPT
                             1 INFEASIBLE (INFES)
                    SUM      9.358
                    MAX      9.358
                    MEAN     9.358
                             0  UNBOUNDED

Cheers
Renger

moatta · October 16, 2019, 12:29pm

Thanks Renger for your reply.

I did not notice the message saying “Model has been proven infeasible” in the other window. I changed =e= to =g= in const3 and I got a reasonable value for the variable OF which is 1056.250. In the past case, I got 9.358 which is impossible because the cheapest MW costs 10 and I need 100 MW. However, I got the same message saying “Model has been proven infeasible”.

So, even =g= gives a wrong answer, right? What do you recommend to get better results?

Renger · October 16, 2019, 3:59pm

Hi
Hard to say: your model is infeasible so it might have no solution with the specifications you give. The constraint 7 seems to be a problem.
You should btw always check the model status in the solve summary:

S O L V E      S U M M A R Y

     MODEL   loadflow            OBJECTIVE  OF
     TYPE    LP                  DIRECTION  MINIMIZE
     SOLVER  CPLEX               FROM LINE  70

**** SOLVER STATUS     1 Normal Completion         
**** MODEL STATUS      4 Infeasible                
**** OBJECTIVE VALUE                9.3584

Cheers
Renger

moatta · November 3, 2019, 9:41am

Thanks Renger for your help. It is working very well now.

Thanks again! I really appreciate it

falfaisa · November 5, 2019, 8:16am

Hello,

I have run your code after disables constraints 7. I got a feasible solution with the OF using =e=. The OF is 1056.250 .

Set
   bus        / 1*3   /
   slack(bus) / 3     /
   Gen        / g1*g3 /;

Alias (bus,node);

Table GenData(Gen,*) 'generating units characteristics'
       b   pmin  pmax
   g1  10  0     65
   g2  11  0     100;
* -----------------------------------------------------

Set GBconect(bus,Gen) 'connectivity index of each generating unit to each bus' / 1.g1, 3.g2 /;

Table BusData(bus,*) 'demands of each bus in MW'
       Pd
   2  100;

Set conex 'bus connectivity matrix' / 1.2, 2.3, 1.3 /;
conex(bus,node)$(conex(node,bus)) = 1;

Table branch(bus,node,*) 'network technical characteristics'
        x     Limit
   1.2  0.2   50
   2.3  0.25  100
   1.3  0.4   100 ;

branch(bus,node,'x')$(branch(bus,node,'x')=0)         =   branch(node,bus,'x');
branch(bus,node,'Limit')$(branch(bus,node,'Limit')=0) =   branch(node,bus,'Limit');
branch(bus,node,'bij')$conex(bus,node)                = 1/branch(bus,node,'x');
*****************************************************

Variable OF, Pij(bus,node), Pg(Gen), delta(bus);
Equation const1, const2, const3, const9, const10, const4, const5, const6,
* const7,
const8;


const1(bus,node)$conex(bus,node)..
   Pij(bus,node) =e= branch(bus,node,'bij')*(delta(bus) - delta(node));

const2(bus)..
   sum(Gen$GBconect(bus,Gen), Pg(Gen)) - BusData(bus,'Pd') =e= sum(node$conex(node,bus), Pij(bus,node));

const3..
   OF =e= sum(Gen,Pg(Gen)*GenData(Gen,'b'));

const9(bus,node)$conex(bus,node)..
         Pij(bus,node)=l=branch(bus,node,'Limit');

const10(bus,node)$conex(bus,node)..
         Pij(bus,node)=g=-branch(bus,node,'Limit');

const4(Gen)..
         Pg(Gen)=l=GenData(Gen,'Pmax');

const5(Gen)..
         Pg(Gen)=g=GenData(Gen,'Pmin');

const6(bus)..
         delta(bus)=l=pi;

*const7(bus)..
*         delta(bus)=g=-pi;

const8(bus)..
         delta('3')=e=0;

Model loadflow / all /;

solve loadflow minimizing OF using lp;
display Pij.l, of.l, pg.l, delta.l;

Good Luck.
Faisal