Getting variable values in PowerModels.jl

Does anyone know how I can get the value of the lifted variables (e.g., _td) in PowerModels.jl? Thanks in advance!

This function may provide a helpful example, PowerModels.jl/obbt.jl at master · lanl-ansi/PowerModels.jl · GitHub

Thanks for the response! In obbt.jl code a variable optimizes first and its optimal value obtains by “getobjectivevalue()” command. I need to minimize generation cost and obtain all lifted variables (e.g., _d, w_fr, etc). Can you please guide me how I can obtain those variables? Is the any specific command that can do this?Thanks in advance!

As of PowerModels v0.8, you can do something like this,

using JuMP; using PowerModels; using Ipopt
pm = build_generic_model("case3.m", QCWRPowerModel, PowerModels.post_opf)
solve_generic_model(pm, IpoptSolver())
getvalue(var(pm, :td))

Thanks for your response. I could obtain every single variable’s vale in the model. However, I need their values to four digits after decimal points. I used the “floor” command but I got an error.

floor(getvalue(var(pm, :td)))

ERROR: MethodError: no method matching floor(::JuMP.JuMPArray{Float64,1,Tuple{Base.KeyIterator{Dict{Tuple{Int64,Int64},Dict{String,Any}}}}})

I will appreciate it if you guide me to solve the issue!

Values of all angle differences are zero for all branches and values of voltage magnitude are 1 for all buses in “nesta_case5_pjm”. I think these values have been rounded before storing in the dictionary. Actual value of these variables (i.e., these values in the optimal solution) are different than those obtained by “getvalue” command. Please let me know if I have done something wrong. I implemented the “getvalue” command like the command in your previous post but I don’t know why the results are not compatible with the optimal solution. Thanks in advance!

If I run,

using JuMP; using PowerModels; using Ipopt
pm = build_generic_model("nesta_case5_pjm.m", QCWRPowerModel, PowerModels.post_opf)
solve_generic_model(pm, IpoptSolver())

getvalue(var(pm, :td))
getvalue(var(pm, :vm))

I have an output like,

julia> getvalue(var(pm, :td))
0_1_td: 1 dimensions:
[(1, 2)] = 0.0737577421787345
[(2, 3)] = 0.0009511415894786759
[(1, 4)] = 0.05449794187473173
[(1, 5)] = -0.013965702498034971
[(3, 4)] = -0.02021094189348145
[(4, 5)] = -0.06846364437276671

julia> getvalue(var(pm, :vm))
0_1_vm: 1 dimensions:
[4] = 1.0948294398423672
[2] = 1.0857389679556273
[3] = 1.0942198705117394
[5] = 1.1
[1] = 1.0999999924705084

For the floor thing you could try something like,

for bp in ids(pm, :buspairs)
    println(floor(getvalue(var(pm, :td)[bp])))
end

Thank you so much for the response! It was useful for me.
I got the same results as you. Voltage angle difference values are appeared by their corresponding (f_bus, t_bus) tuples. Is there any way to put voltage angle difference values and their corresponding tuples (i.e., f_bus, t_bus) together in a matrix and sort rows of the matrix based on voltage angle difference values. I need that to determine branches with large angle differences. I will appreciate your response!

This is probably how I would compute that (tested in Julia v0.6),

case = PowerModels.parse_file("nesta_case5_pjm.m")
result = run_opf(case, ACPPowerModel, IpoptSolver())
PowerModels.update_data(case, result["solution"])

bus_va = Dict( parse(Int,i) => bus["va"] for (i,bus) in case["bus"])
branch_vad = Dict(i => bus_va[branch["f_bus"]] - bus_va[branch["t_bus"]] for (i,branch) in case["branch"])

branch_order = sort(collect(branch_vad), by=x -> -abs(x[2]))

The advantage of not digging into the model variables is that it can work with a wider variety of power formulations.

Thanks for the response! Just one more question, what is the best way to convert bus_va or bra nch_order elements to Float or Int. I used the following command to convert the keys in bus_va into Float but I couldn’t convert values in bus_va into Float. Can you please let me know how I can do that? Note that I have to use “imap” command because I need orders to be preserved. Thanks in advance!

collect(imap(first,bus_va))

Does this resolve your issue?

bus_va = Dict( parse(Int,i) => bus["va"] for (i,bus) in case["bus"])
branch_vad = Dict( parse(Int,i) => bus_va[branch["f_bus"]] - bus_va[branch["t_bus"]] for (i,branch) in case["branch"])

branch_order = sort(collect(branch_vad), by=x -> -abs(x[2]))

No, it did resolve the issue to some extent. However, output for these command are still dictionary which I cannot turn them into matrices easily. For instance following is “branch_order”

6=>-0.0626634
1=>0.0617572
2=>0.0489352
3=>-0.0137282
5=>-0.00976905
4=>-0.00305291

But I need that to be in matrix form

[ 6 , -0.0626634
1 , 0.0617572
2 , 0.0489352
3 ,-0.0137282
5 ,-0.00976905
4 ,-0.00305291]

I will appreciate it if you guide me to turn them to matrices. Thanks in advance!

I used following commands to change the dictionary first into a tuple and then change the tuple into array.

a = [(k,v) for (k,v) in branch_order]
t = [y[i] for y in a, i in 1:2]

It works for me but I surfed the net and people were pessimistic about performance of these commands. Do you think these commands make the code sluggish? Or do have a better suggestion? Thanks in advance!

You can just go

[key_value[i] for key_value in branch_order, i in 1:2]

Thanks for the response!
Following is the results that I have obtained for a variable.

[(1, 2)] = 0.0737577421787345
[(2, 3)] = 0.0009511415894786759
[(1, 4)] = 0.05449794187473173
[(1, 5)] = -0.013965702498034971
[(3, 4)] = -0.02021094189348145
[(4, 5)] = -0.06846364437276671

Can you please guide me how I can extract just the values? Thanks in advance!

0.0737577421787345
0.0009511415894786759
0.05449794187473173
-0.013965702498034971
-0.02021094189348145
-0.06846364437276671

I would not worry too much about performance until you run into runtime issues. Most power network datasets are relatively small.

Thanks for the response! When I executed “getvalue(var(pm, :td))” command I get the following results.
[(1, 2)] = 0.0726892490132161
[(2, 3)] = 0.0007423729745521509
[(1, 4)] = 0.05311970563797517
[(1, 5)] = -0.013538436007847555
[(3, 4)] = -0.020311916349793093
[(4, 5)] = -0.06665814164582272

but when I run " branch_vad = Dict( parse(Int,i) => bus_va[branch[“f_bus”]] - bus_va[branch[“t_bus”]] for (i,branch) in case[“branch”])" I get the following.

4 => 0.000742373
2 => 0.0531197
3 => -0.0135384
5 => -0.0203119
6 => -0.0666581
1 => 0.0726892

As you can see branch indices in the outputs are not consistent. I do want to compare voltage angle difference variable (td) with its corresponding (va(f_bus)-va(t_bus)) value. I will appreciate it if you let me know about a way to rearrange one of these outputs to have a same branch index with another one. Thanks in advance and I am looking forward to hearing from you.

There is a many-to-one relationship between the indices of :td and the branch ids (e.g. consider parallel branches).

My recommendation would be to build your code around the branch ids used in the later example. If you also want easy access to the bus ids, give this a try,

branch_vad = Dict( (parse(Int,i), branch[“f_bus”], branch[“t_bus”])  => bus_va[branch[“f_bus”]] - bus_va[branch[“t_bus”]] for (i,branch) in case[“branch”])

Thanks for the instructive recommendation. I took your word and put the piece of code that returns voltage angle difference values, in a for loop in order to have voltage angle difference for all branches including parallel ones. The issue is resolved and now voltage angle difference array obtain by “td” and (“va(f_bus)-va(t_bus)”) have a same length (i.e., the length of branch matrix).

aa = getvalue(var(pm, :td))
for (l,branch) in case[“branch”]
i = branch[“f_bus”]
j = branch[“t_bus”]
bpp= Dict( (parse(Int,i), branch[“f_bus”] , branch[“t_bus”]) => aa[(branch[“f_bus”],branch[“t_bus”])] for (i,branch) in case[“branch”])
end

Thank you so much once again!

I have an array (w) including binaries that determines weather or not a specific action should taken on each branch. In other words if w[i] =1 then the action would applicable on ith branch. For instance, for “nesta_case5_pjm” I have w=[0;1;1;0;0;0]. I want to convert this array to a dictionary in which values in s can be callable by “bp”. I tried to do that by following command

bpp= Dict( (branch["f_bus"] , branch["t_bus"]) => w[i] for (i,branch) in case["branch"])    

but I got the following error. Can you please guide me on doing this? Thanks in advance!

ERROR: ArgumentError: invalid index: 4