Coding up basic power flow equation

mbah · September 8, 2021, 8:18pm

Hello,

I am relatively new to Julia and power systems optimization. At the moment, I am trying to understand power flows from scratch and this involves coding up/defining branch power flows for a simple
three bus network.

3node example

To do this, I have resorted to a crude approach of first defining the line equations and then attaching the equation to an equality constraint. I also want to limit the power flows to some value. The approach isn’t working as I’m getting error messages.


using JuMP, Clp, DataFrames, CSV

basic_model = Model(Clp.Optimizer)

@variable(basic_model, GEN[T, C] >=0) #Conventional plant generation
@variable(basic_model, FLOW[T, L] ) #power flow

@constraint(basic_model, 
  CapConstraint[t in T, c in C],
  GEN[t,c] <= data[:plants][:g_max][c] ) # Capacity constraint

## issue starts here

flow1 = sum(1/3*GEN[1:3,"basel"] - 1/3*GEN[1:3,"bern"])
flow2 = sum(2/3*GEN[1:3,"basel"] + 1/3*GEN[1:3,"bern"])
flow3 = sum(1/3*GEN[1:3,"basel"] + 2/3*GEN[1:3,"bern"])

@constraint(basic_model, 
    LineFlowA[t in T, l in L],
    FLOW[t,l] == flow1 )

@constraint(basic_model, 
    LineFlowB[t in T, l in L],
    FLOW[t,l] == flow2 )

@constraint(basic_model, 
    LineFlowC[t in T, l in L],
    FLOW[t,l] == flow3 )

@constraint(basic_model, 
    FlowMax_Pos[t in T, l in L],
    FLOW[t,l] <= data[:lines][:f_max][l]
    ) #positve line flow  constraint

@constraint(basic_model, 
    FlowMax_Neg[t in T, l in L],
    FLOW[t,l] >= -data[:lines][:f_max][l]
    ) #negative line flow constraint

Sample error message when trying to code the flow (flow1 - flow3)


ERROR: MethodError: no method matching promote_shape(::Tuple{Vector{Int64}}, ::Tuple{Vector{Int64}})

A simple and viable approach to define the line flow constraints would be greatly appreciated.

Thanks

ccoffrin · September 8, 2021, 8:33pm

Hi @mbah, it’s nice to hear you are giving Julia and JuMP a try for the power flow problem.

I can recommend these reference implementations for solving DC and AC optimal power flow,

https://github.com/lanl-ansi/PowerModelsAnnex.jl/tree/master/src/model

Adapting these examples into a power flow formulation is a good exercise.

You may also find the power flow solvers in PowerModels.jl to be a useful reference for other approaches to modeling the power flow problem in Julia.

https://lanl-ansi.github.io/PowerModels.jl/stable/power-flow/

mbah · September 8, 2021, 8:53pm

Thanks for the super quick response @ccoffrin. I will follow up on your suggestion.

odow · September 9, 2021, 12:50am

Hi there! Since this is your first post, welcome.

To add on to what @ccoffrin said: I’m not sure what the problem is because I can’t run your code. We might be able to provide more insight if you can provide the (simplified) data, or if you can provide the full error message. This post has some helpful tips on writing a good question: Please read: make it easier to help you

mbah · September 9, 2021, 9:49pm

Thanks for the response @odow.
Since, I’m unable to upload the csv files here, I have made the simplified data available on my github page: https://github.com/m-bah/test

A slightly more expanded version of my coding exercise;

using JuMP, Clp, DataFrames, CSV
using Base: Symbol

function df_to_dict_with_id(df::DataFrame)
  Dict(col[1] => Dict(zip(df[!, 1], col[2])) for col in pairs(eachcol(df)))
end

#Read in data

data_dir = pwd()*"/test/"
data = Dict{Symbol,Any}()
set = Dict{Symbol,Any}() 
readDataDFs = Dict()
inputDataItems = ["plants", "demand", "nodes","lines"]

for inputDataItem in inputDataItems #loop starts~ for each data file inputDataItems
  readDataDFs[Symbol(inputDataItem)] = DataFrame(CSV.File(data_dir * inputDataItem * ".csv"), copycols=true) #read in all the CSV files into the dictionary readDataDFs
  data[Symbol(inputDataItem)] = df_to_dict_with_id(readDataDFs[Symbol(inputDataItem)]) 
  set[Symbol(inputDataItem)] = readDataDFs[Symbol(inputDataItem)].index 
end

P = set[:plants] #index of all plants
C = filter(r -> any(occursin.(["conventional"], r.plant_type)), readDataDFs[:plants]).index #index of conventional plants
N = set[:nodes] #Index of nodes  
T = set[:demand] #Time stamp
L = set[:lines] #Index of power lines

basic_model = Model(Clp.Optimizer) 

#Define variables
@variable(basic_model, GEN[T, C] >=0) #Conventional plant generation
@variable(basic_model, FLOW[T, L] ) #power flow

@objective(basic_model, 
  Min, 
  sum(sum(GEN[t, p]*data[:plants][:cost][p] for p in C) for t in T)) 

@constraint(basic_model, 
  CapConstraint[t in T, c in C],
  GEN[t,c] <= data[:plants][:g_max][c] ) # Capacity constraint

.... #steps skipped

## issue starts here

flow1 = sum(1/3*GEN[1:3,"basel"] - 1/3*GEN[1:3,"bern"])
flow2 = sum(2/3*GEN[1:3,"basel"] + 1/3*GEN[1:3,"bern"])
flow3 = sum(1/3*GEN[1:3,"basel"] + 2/3*GEN[1:3,"bern"])

@constraint(basic_model, 
    LineFlowA[t in T, l in L],
    FLOW[t,l] == flow1 )

@constraint(basic_model, 
    LineFlowB[t in T, l in L],
    FLOW[t,l] == flow2 )

@constraint(basic_model, 
    LineFlowC[t in T, l in L],
    FLOW[t,l] == flow3 )

@constraint(basic_model, 
    FlowMax_Pos[t in T, l in L],
    FLOW[t,l] <= data[:lines][:f_max][l]
    ) #positve line flow  constraint

@constraint(basic_model, 
    FlowMax_Neg[t in T, l in L],
    FLOW[t,l] >= -data[:lines][:f_max][l]
    ) #negative line flow constraint

When I try to compute flow1 for instance, i get the following error message.

julia> flow1 = sum(1/3*GEN[1:3,"basel"] - 1/3*GEN[1:3,"bern"])
ERROR: MethodError: no method matching promote_shape(::Tuple{Vector{Int64}}, ::Tuple{Vector{Int64}})
Stacktrace:
 [1] promote_shape(a::JuMP.Containers.DenseAxisArray{AffExpr, 1, Tuple{Vector{Int64}}, Tuple{JuMP.Containers._AxisLookup{Dict{Int64, Int64}}}}, b::JuMP.Containers.DenseAxisArray{AffExpr, 1, Tuple{Vector{Int64}}, Tuple{JuMP.Containers._AxisLookup{Dict{Int64, Int64}}}})
   @ Base .\indices.jl:169
 [2] -(A::JuMP.Containers.DenseAxisArray{AffExpr, 1, Tuple{Vector{Int64}}, Tuple{JuMP.Containers._AxisLookup{Dict{Int64, Int64}}}}, B::JuMP.Containers.DenseAxisArray{AffExpr, 1, Tuple{Vector{Int64}}, Tuple{JuMP.Containers._AxisLookup{Dict{Int64, Int64}}}})
   @ Base .\arraymath.jl:38
 [3] top-level scope
   @ REPL[26]:1

In short, my main aim is to manually code up equality constraints for individual line flows (i.e. LineFlowA - LineFlowC).

Apologies for the long and messy response

Thanks

odow · September 9, 2021, 10:49pm

Try:

flow1 = sum(1/3*GEN[1:3,"basel"] .- 1/3*GEN[1:3,"bern"])

You want element-wise subtraction, so you need .- instead of -.

mbah · September 15, 2021, 9:55am

Thanks a lot @odow for the suggestion. It worked out!
Below is a snippet of the solution I used to hard-code individual line constraints.

flow1 = 1/3*GEN[1:3,"basel"] .- 1/3*GEN[1:3,"bern"]
flow2 = 2/3*GEN[1:3,"basel"] .+ 1/3*GEN[1:3,"bern"]
flow3 = 1/3*GEN[1:3,"basel"] .+ 2/3*GEN[1:3,"bern"]


@constraint(basic_model,
    LineFlowA[t in T],
      FLOW[t,"BS-BE"] == flow1[t]) 

@constraint(basic_model, 
    LineFlowB[t in T],
    FLOW[t, "BS-ZH"] == flow2[t])

@constraint(basic_model, 
    LineFlowC[t in T],
    FLOW[t,"BE-ZH"] == flow3[t])

Unfortunately, the final solution is infeasible.

julia> optimize!(basic_model)
Coin0507I Presolve determined that the problem was infeasible with tolerance of 1e-08
Clp3003W Analysis indicates model infeasible or unbounded
Clp0006I 0  Obj 0 Primal inf 200 (3)
Clp0006I 3  Obj 0 Primal inf 300 (4)
Clp0006I 5  Obj 3000 Primal inf 500 (4)
Clp0001I Primal infeasible - objective value 3000
Clp0032I PrimalInfeasible objective 3000 - 5 iterations time 0.022

I changed the solver from Clp to Ipopt and I get the same results, although with a bit more information as to why the solution is infeasible

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.13.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:       51
Number of nonzeros in inequality constraint Jacobian.:        6
Number of nonzeros in Lagrangian Hessian.............:        0

Exception of type: TOO_FEW_DOF in file "IpIpoptApplication.cpp" at line 938:
 Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!

EXIT: Problem has too few degrees of freedom.

Any ideas/solutions to this problem would be greatly appreciated.

Thanks