# Newton solution from PowerModels.jl

**URL:** <https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356>\
**Category:** Optimization (Mathematical)\
**Created:** [July 9, 2021, 2:08pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356 "2021-07-09T14:08:25Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mariana](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mariana/32/15064_2.png) [@Mariana](https://discourse.julialang.org/u/Mariana)\
**Post date:** [July 9, 2021, 2:08pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/1 "2021-07-09T14:08:25Z")

</div>

Hello.

I’m using PowerModels.jl to run some power flow studies and I’m interested in comparing the results obtained from the JuMP solution to the Newton solution. In most cases, those results are the same; however, when I make modifications on the power system (e.g. deleting a transmission line) and run Newton’s method, the flow on one branch of my system results in a much higher value when compared to JuMP’s solution. I’m using version 0.18.1 of PowerModels and my system has 10 buses.

This is the code I’m using to run this power flow study:

```julia
using PowerModels
using Ipopt

network_data = parse_file("/teste_novo.raw")
print_summary(network_data)

# ========================= Changes in the system =========================

# Open line 8-9 (branch 6):
delete!(network_data["branch"], "6")

# ========================= Newton solution =========================

result = compute_ac_pf(network_data)

# Check that the solver converged:
update_data!(network_data, result["solution"])

flows = calc_branch_flow_ac(network_data)
update_data!(network_data, flows)
print_summary(network_data)

```

Am I doing something wrong in this code for Newton’s method?

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [July 9, 2021, 2:33pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/2 "2021-07-09T14:33:53Z")

</div>

Hi @Mariana, the work flow you are showing here looks reasonable. So I it is not clear what might be the problem to me. One question that does come to mind, when you speak of the “JuMP solution”, I presume you are referring to the output of a call to `run_ac_pf`?

---

<div class="post-metadata">

**Author:** ![Mariana](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mariana/32/15064_2.png) [@Mariana](https://discourse.julialang.org/u/Mariana)\
**Post date:** [July 9, 2021, 2:44pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/3 "2021-07-09T14:44:39Z")

</div>

Yes! In fact, I’m using the `run_pf` function as:

```julia
run_pf(network_data, ACPPowerModel, Ipopt.Optimizer)

```

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [July 9, 2021, 4:43pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/4 "2021-07-09T16:43:22Z")

</div>

Ok. It is not obvious to me why you would see different results in these two cases. It is possible the two approaches have converged to different feasible flow solutions, but this is very rare in my experience. This kind of workflow can be used to test the feasibility of a given power flow solution,

> <https://github.com/lanl-ansi/PowerModels.jl/blob/master/test/data.jl#L787>

If you find there are very large power imbalances this might be a bug in PowerModels, but I would need a detailed example to confirm and fix it, if that is the case.

---

<div class="post-metadata">

**Author:** ![Mariana](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mariana/32/15064_2.png) [@Mariana](https://discourse.julialang.org/u/Mariana)\
**Post date:** [July 16, 2021, 9:56am UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/5 "2021-07-16T09:56:52Z")

</div>

@ccoffrin sorry for the late response! I followed the workflow you posted and, in order to eliminate branch 6 of my system, I ran the following code:

```julia
@testset "test 3" begin
    data = PowerModels.parse_file("/teste_novo.raw")
    data["branch"]["6"]["br_status"] = 0
    result = compute_ac_pf(data)
    #result = run_pf(data, ACPPowerModel, Ipopt.Optimizer)
    PowerModels.update_data!(data, result["solution"])

    flows = PowerModels.calc_branch_flow_ac(data)
    PowerModels.update_data!(data, flows)

    balance = PowerModels.calc_power_balance(data)

    for (i,bus) in balance["bus"]
        @test isapprox(bus["p_delta"], 0.0; atol=1e-6)
        @test isapprox(bus["q_delta"], 0.0; atol=1e-6)
    end
end

```

All the 20 tests pass when I use the `run_pf` function; however, when I use the `compute_ac_pf` function, 16 tests pass and 4 fail.

Am I doing something wrong?

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [July 16, 2021, 2:44pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/6 "2021-07-16T14:44:41Z")

</div>

I would start by checking the status value in `result` to see if `compute_ac_pf` belives it has converged or not. You can also pass NLsolve keyword arguments to `compute_ac_pf` to control the options of that solver.

> **[GitHub - JuliaNLSolvers/NLsolve.jl: Julia solvers for systems of nonlinear...](https://github.com/JuliaNLSolvers/NLsolve.jl)**
>
> Julia solvers for systems of nonlinear equations and mixed complementarity problems - GitHub - JuliaNLSolvers/NLsolve.jl: Julia solvers for systems of nonlinear equations and mixed complementarity ...

I will also point you to this documentation about the different power flow models that are available in PowerModels,

[https://lanl-ansi.github.io/PowerModels.jl/stable/power-flow/](https://lanl-ansi.github.io/PowerModels.jl/stable/power-flow/)

`compute_ac_pf` is not expected to be as reliable as `run_ac_pf` using a solver like Ipopt.

---

<div class="post-metadata">

**Author:** ![Mariana](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mariana/32/15064_2.png) [@Mariana](https://discourse.julialang.org/u/Mariana)\
**Post date:** [July 22, 2021, 5:22pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/7 "2021-07-22T17:22:33Z")

</div>

@ccoffrin thank you very much for your answers! I tried the workflow you suggested to check the feasibility of the solution with different power systems and all tests passed, so I guess I was using an ill-conditioned system. Just to confirm: the `PowerModels.calc_power_balance` function works by checking if Kirchhoff’s Law is satisfied in all the buses of the system, right?

---

<div class="post-metadata">

**Author:** ![Mariana](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mariana/32/15064_2.png) [@Mariana](https://discourse.julialang.org/u/Mariana)\
**Post date:** [July 22, 2021, 5:26pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/8 "2021-07-22T17:26:01Z")

</div>

Another question: could you explain, or point to some documentation, how PowerModels.jl solve power flows with multiple reference buses? My system has 9 buses, and 2 reference buses. When I run `run_ac_pf` and `compute_ac_pf` I get the same results, but they are different than the ones I obtained from another power flow software.

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [July 22, 2021, 8:02pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/9 "2021-07-22T20:02:37Z")

</div>

> [@Mariana](#):
>
> Just to confirm: the `PowerModels.calc_power_balance` function works by checking if Kirchhoff’s Law is satisfied in all the buses of the system, right?

Yes that is correct. By default it does not calculate the branch flow values so before calling `calc_power_balance` you should make sure the branch flows have been computed using the latest version of the voltage profile.

> [@Mariana](#):
>
> how PowerModels.jl solve power flows with multiple reference buses?

I am not aware of how other software treats multiple reference buses, so I am not surprised there may be some differences here.

In the case of `compute_ac_pf` I think this should produce an error or warning message. If not, please post an issue.

For `run_ac_pf` PowerModels will build one JuMP model that has two PowerFlow models inside of it. As long as they are mathematically well formed it should converge. Typically “mathematically well formed” AC Power Flow requires that there is exactly one reference bus per connected component in the network. The function `correct_bus_types!` can be used to check and prepare a given network data for Power Flow solving.

> <https://github.com/lanl-ansi/PowerModels.jl/blob/master/src/core/data.jl#L1687>

---

<div class="post-metadata">

**Author:** ![iagoschavarry](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iagoschavarry/32/27460_2.png) [@iagoschavarry](https://discourse.julialang.org/u/iagoschavarry)\
**Post date:** [July 23, 2021, 2:21pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/10 "2021-07-23T14:21:55Z")

</div>

Hi @ccoffrin!

> [@ccoffrin](#):
>
> In the case of `compute_ac_pf` I think this should produce an error or warning message. If not, please post an issue.

I think since 0.3.4 the code is generic for any number of reference buses when running AC Power Flow. If more than one reference is given, the angle of the multiple references buses are set to zero. I also don’t know how other softwares treat this feature.

Currently, `correct_bus_types!` only checks if there is an active generator in each given reference bus, but it does not limit the number of ref buses to one.

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [July 23, 2021, 4:00pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/11 "2021-07-23T16:00:28Z")

</div>

> [@iagoschavarry](#):
>
> Currently, `correct_bus_types!` only checks if there is an active generator in each given reference bus, but it does not limit the number of ref buses to one.

This is correct, thanks for that clarification @iagoschavarry!

> [@iagoschavarry](#):
>
> I think since 0.3.4 the code is generic for any number of reference buses when running AC Power Flow.

In the case of `run_ac_pf` you are correct. `compute_ac_pf` is a relatively new feature that solves an AC power flow without building a JuMP model, so its functionalities are more limited. You can see in this test script that multiple slack buses are currently not supported.

> <https://github.com/lanl-ansi/PowerModels.jl/blob/master/test/pf-native.jl#L143>

---

<div class="post-metadata">

**Author:** ![iagoschavarry](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iagoschavarry/32/27460_2.png) [@iagoschavarry](https://discourse.julialang.org/u/iagoschavarry)\
**Post date:** [July 23, 2021, 6:08pm UTC](https://discourse.julialang.org/t/newton-solution-from-powermodels-jl/64356/12 "2021-07-23T18:08:00Z")

</div>

> [@ccoffrin](#):
>
> In the case of `run_ac_pf` you are correct. `compute_ac_pf` is a relatively new feature that solves an AC power flow without building a JuMP model, so its functionalities are more limited. You can see in this test script that multiple slack buses are currently not supported.

Right, @ccoffrin! Thanks for the response.

I tried to follow the begin-to-end code of `compute_ac_pf`.

In 0.18.2 it seems that this function is also generic to any numbers of reference buses. Although, this appers to be a recent feature, added in the 0.18.0 update.

Before that, 0.17.4 for example, `compute_ac_pf` calls `calc_admittance_matrix` which calls `reference_bus`. The last function checks if there is more than one reference bus in the network data.  
Since 0.18.0, `reference_bus` is no longer called in `calc_admittance_matrix`, therefore no warnings or errors are given when multiple references buses are used in `compute_ac_pf`.

I guess this wasn’t a intended feature then? In our studies, calling `compute_ac_pf` and `run_ac_pf` with mutiple reference buses in a small network data (10 buses) resulted approximately in the same outcome.