Python - Help with constraints

1,053 views
Skip to first unread message

Juan Zhang

unread,
Dec 4, 2017, 12:21:31 AM12/4/17
to Gurobi Optimization
Can anyone help me figure out what is wrong with my syntax on the constraints?

Here is the python Code:
from gurobipy import *

Demand = [0, 300, 300, 300, 0]

GenToBusMatrix = [
    [1, 0, 0, 0, 0],
    [1, 0, 0, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 0, 1]
]

LineToBusMatrix = [
    [1, -1, 0, 0, 0],
    [1, 0, 0, -1, 0],
    [0, 1, -1, 0, 0],
    [0, 0, 1, -1, 0],
    [0, 0, 0, 1, -1],
    [1, 0, 0, 0, -1]
]

LineReactance = [100/0.0291, 100/0.0304, 100/0.0108, 100/0.0297, 100/0.297, 100/0.0064]

LinePowerLimit = [1000, 1000, 1000, 1000, 240, 1000]

def SolveEDwithPowerLines (NGenUnits, PMin, PMax, a, b ,c, NBuses, NLines, GeneratorToBusMatrix, LineToBusMatrix1, LineReactance, LinePowerLimit, Demand):

    #####################################################################################################################################
    ################################################        INITIALIZATION       ########################################################
    #####################################################################################################################################
    units = [0] * NGenUnits
    lines = [0] * NLines
    buses = [0] * NBuses

    Generation_Data = {}
    GenToBusMatrix = {}
    LineToBusMatrix = {}
    LineDataMatrix= {}

    for i in range(0, NBuses):
        buses [i] = i
    for i in range(0, NLines):
        lines [i] = i
        LineDataMatrix [i, 'x'] = LineReactance[i]
        LineDataMatrix [i, 'PLmax'] = LinePowerLimit [i]


    for j in range(0, NGenUnits):
        units[j] = j
        Generation_Data [j, 'a'] = a[j]
        Generation_Data [j, 'b'] = b[j]
        Generation_Data [j, 'c'] = c[j]
        Generation_Data [j, 'PMin'] = PMin[j]
        Generation_Data [j, 'PMax'] = PMax[j]

    for i in range(0,NLines):
        for j in range(0,NBuses):
            LineToBusMatrix [i, j] = LineToBusMatrix1[i][j]

    for i in range(0,NBuses):
        for j in range(0,NGenUnits):
            GenToBusMatrix[i, j] = GeneratorToBusMatrix[i][j]

    #####################################################################################################################################
    ############################################        MODEL DEFINITION       ##########################################################
    #####################################################################################################################################

    m = Model("ED_withPowerLines")

    #####################################################################################################################################
    ################################################        VARIABLES       #############################################################
    #####################################################################################################################################

    P =                 m.addVars(units, name="Power_Generation")
    PowerFlow_Lines =   m.addVars(lines, name="Power on the lines")
    theta =             m.addVars(buses, name="Bus Angles")
    theta[0] = 0
    Total_Cost = m.addVar(vtype=GRB.BINARY, name="Total_Cost")

    m.update()

    #####################################################################################################################################
    ############################################        OBJECTIVE FUNCTION       ########################################################
    #####################################################################################################################################

    obj = quicksum(
        Generation_Data[unit, 'a']
        + Generation_Data[unit, 'b'] * P[unit]
        + Generation_Data[unit, 'c'] * P[unit] * P[unit]
        for unit in units
    )

    m.setObjective(obj, GRB.MINIMIZE)

    #####################################################################################################################################
    ################################################        CONSTRAINTS       ###########################################################
    #####################################################################################################################################

    m.addConstr(
        (PowerFlow_Lines[line]  for line in lines) ==
        ((quicksum(theta[bus] * LineToBusMatrix[line, bus] * LineDataMatrix[line,'x']) for line in lines for bus in buses))
    )
    m.addConstr((P[unit] for unit in units) <= (Generation_Data[unit, 'PMax'] for unit in units), name="PMax")
    m.addConstr((P[unit] for unit in units) >= (Generation_Data[unit, 'PMin'] for unit in units), name="PMin")
    m.addConstr(
        (quicksum((P[unit] * GenToBusMatrix[bus, unit]) for bus in buses for unit in units) - (Demand[bus] for bus in buses)) ==
        (quicksum(LineToBusMatrix[line, bus]) for line in lines for bus in buses)
    )

    m.addConstr((PowerFlow_Lines[line] for line in lines) <= (LineDataMatrix[line, 'PLmax'] for line in lines))
    m.addConstr((PowerFlow_Lines[line] for line in lines) >= (-LineDataMatrix[line, 'PLmax'] for line in lines))

    #####################################################################################################################################
    ##############################################        RUN OPTIMIZER       ###########################################################
    #####################################################################################################################################

    m.optimize()

    #####################################################################################################################################
    #############################################        RETRIEVE RESULTS       #########################################################
    #####################################################################################################################################

    for v in m.getVars():
        if v.X != 0:
            print("%s %f" % (v.Varname, v.X))

#########################################################################################################################################
#####################################################        END       ##################################################################
#########################################################################################################################################

SolveEDwithPowerLines(
    5, #Number of Generation Units
    [0, 0 ,0 ,0, 0], #Pmin
    [110, 100, 520, 200, 600], #PMax
    [0, 0, 0, 0, 0], #a
    [14, 15, 30, 29, 10], #b
    [0, 0, 0, 0, 0], #c
    5, # Number of Buses
    6, # #Number of lines
    GenToBusMatrix,
    LineToBusMatrix,
    LineReactance,
    LinePowerLimit,
    Demand
    )

I have solved this in GAMS, here is the code:

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

Sets
NG      set of generators /i1*i5/
NL      set of lines  /l1*l6/
NB      set of buses /b1*b5/
;


table UnitData(NG,*)     Generation units data
      Pmin    Pmax      A       B         C
i1    0       110       0       14        0
i2    0       100       0       15        0
i3    0       520       0       30        0
i4    0       200       0       29        0
i5    0       600       0       10        0
;


Table  UnitToBus(NG,NB)  Mapping from generators to buses
      b1    b2    b3     b4      b5
i1    1     0     0      0       0
i2    1     0     0      0       0
i3    0     0     1      0       0
i4    0     0     0      1       0
i5    0     0     0      0       1
;

Table  LineToBus(NL,NB)  Mapping from lines to buses
      b1    b2    b3     b4      b5
l1    1     -1     0     0       0
l2    1     0      0     -1      0
l3    0     1     -1     0       0
l4    0     0      1     -1      0
l5    0     0      0     1       -1
l6    1     0      0     0       -1
;

Table LineData(NL,*)     Line data
      x         PLmax
l1    0.0291    1000
l2    0.0304    1000
l3    0.0108    1000
l4    0.0297    1000
l5    0.0297    240
l6    0.0064    1000
;



parameter Demand(NB)    Demand of bus b at time t
/
b1  0
b2  300
b3  300
b4  300
b5  0
/
;


Positive Variables
P(NG)        Real power generation of generation unit i at time t
;
Variables
Z
Theta(NB)       Voltage angle of bus b
PL(NL)
;
Theta.up(NB)=3.14;
Theta.lo(NB)=-3.14;
Theta.fx('b1')=0;


Equations
Obj               Objective function
Bal(NB)        Demand-Supply balance
Line_Flow(NL)  Line power flow
Max_PL1(NL)    Maximum line power flow 1
Max_PL2(NL)    Maximum line power flow 2
Max_P(NG)      Maximum generation
Min_P(NG)      Minimum generation
;



Obj..                            Z=e=sum(NG, UnitData(NG,'A') + UnitData(NG,'B')*P(NG) + UnitData(NG,'C')*P(NG)**2);

Line_Flow(NL)..                  PL(NL)=e=100*sum(NB,Theta(NB)*LineToBus(NL,NB))/LineData(NL,'x');

Bal(NB)..                        sum(NG,UnitToBus(NG,NB)*P(NG))-Demand(NB)=e=sum(NL,LineToBus(NL,NB)*PL(NL));

Max_PL1(NL)..                    PL(NL)=l=LineData(NL,'PLmax');

Max_PL2(NL)..                    PL(NL)=g=-LineData(NL,'PLmax');

Max_P(NG)..                      P(NG)=l=UnitData(NG,'Pmax');

Min_P(NG)..                      P(NG)=g=UnitData(NG,'Pmin');

model ED /all/;

solve ED minimizing Z using NLP;

Tobias Achterberg

unread,
Dec 4, 2017, 7:10:19 AM12/4/17
to gur...@googlegroups.com
Let's take a look at your last constraint:

m.addConstr((PowerFlow_Lines[line] for line in lines) >= (-LineDataMatrix[line,
'PLmax'] for line in lines))

What should this constraint do? PowerFlow_Lines is a dictionary of Gurobi variables,
LineDataMatrix is a two-dimensional array. Do you want to get a lower bound for each
individual variable? Then you need to add multiple constraints like this:

for line in lines:
m.addConstr(PowerFlow_Lines[line] >= -LineDataMatrix[line, 'PLmax'])

or alternatively use the addConstrs() method to add them in a single shot:

m.addConstrs((PowerFlow_Lines[line] >= -LineDataMatrix[line, 'PLmax']) for line in lines)

I hope this helps to understand also the remaining issues in your code.


Regards,

Tobias

Juan Zhang

unread,
Dec 27, 2017, 2:09:28 AM12/27/17
to Gurobi Optimization
Thank you for your help. 

You are correct, I am trying to set upper and lower bounds for each variable. 

by hanging my syntax to your suggestion I got the following error:

m.addConstr(((P[unit] <= Generation_Data[unit, 'PMax']) for unit in units), name="PMax")
  File "model.pxi", line 2718, in gurobipy.Model.addConstr (../../src/python/gurobipy.c:81435)
TypeError: unsupported operand type(s) for -: 'generator' and 'NoneType'

Do you have any suggestions?

Michael Winkler

unread,
Dec 27, 2017, 5:44:04 AM12/27/17
to Gurobi Optimization
I believe either P[unit] or Generation_Data[unit, 'PMax'] in 'None', so one of them is not set correctly.

Best,
Michael

Tobias Achterberg

unread,
Dec 27, 2017, 8:23:34 AM12/27/17
to gur...@googlegroups.com
Sorry, I don't know. Which exact Gurobi version are you using?

Tobias

Juan Zhang

unread,
Dec 27, 2017, 7:02:52 PM12/27/17
to Gurobi Optimization
I am on gurobipy version 7.5.1, and I am using python 2.7

Wang Minke

unread,
Dec 28, 2017, 8:51:37 AM12/28/17
to Gurobi Optimization
what's your coding editor? Thanks!

Juan Zhang

unread,
Dec 28, 2017, 11:29:55 AM12/28/17
to Gurobi Optimization
I use Pycharm
Reply all
Reply to author
Forward
0 new messages