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

[moatta]

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 / 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]

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

Code: Select all

 and 

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

Cheers
Renger

[moatta]

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

Code: Select all

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]

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

Code: Select all

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

Code: Select all

---- 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]

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]

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:

Code: Select all

               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]

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

Thanks again! I really appreciate it

[falfaisa]

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 .

Code: Select all

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