Maybe Manopt?

[kellertuer]

For now this is just a first sketch of an MWE, because neither Optim nor the MWEs here gave me enough understanding in what they actually solve. None of the variables is speaking, none of the functions is documented – so after about an hour I gave up trying to reverse-engineer what even the 2D example should do.

So this MWE for now just illustrates two quasi Newton Calls and has a lot of comments on what ALM might do.

If we would get an understanding of the constraints, it would be best to implement each (that is each real valued) constraint as its own function, the same for their gradients and pass these to g= and grad_g= of the augmented Lagrangian, cf the docs at

https://manoptjl.org/stable/solvers/augmented_Lagrangian_method/

Things to keep in mind: This will probably just take last place, for a very simple reason – constraint optimisation is – as far as I am aware of – not that old on manifolds yet, see the paper where the ALM is from https://link.springer.com/article/10.1007/s00245-019-09564-3 its from 2020. So this solver does have manifolds in mind in its generality and compared to IPNweton or such (faster and more tailored towards Euclidean problems), this will probably just be super slow.

Also box constraints or such are not socially treated – for actually the same reason, that all this is relatively new stuff on Manifolds. So for $a \leq x \ leq b$ we would actually need two constraints to be implemented. There is also no tools (unless you use some AD from somewhere else) to give you the gradients. Even not for linear constraints, since on Manifolds – linearity does not exist (unless you are Euclidean)

If this is too rough in an idea to add Manopt.jl I can also understand if we just close this.

[odow]

Let me take a look

[kellertuer]

@ccoffrin asked me to try it, I am however – as one might see – not convinced this effort will be useful, the conclusion is probably, as long as the problem is not defined on a manifold (for example the sphere) – there is far better packages around for the Euclidean cases considered here.

[odow]

Here's where I got to. But I waited a while and it didn't even finish the case3 example.

I guess the need to explicit AD (we can do better than ForwardDiff) and the fact we need variable bounds as constraints. But yeah from looking around, I don't know if this type of problem is what Manopt was designed for.

import PowerModels
using Manopt, Manifolds

function solve_opf(file_name)
    time_data_start = time()

    data = PowerModels.parse_file(file_name)
    PowerModels.standardize_cost_terms!(data, order=2)
    PowerModels.calc_thermal_limits!(data)
    ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

    data_load_time = time() - time_data_start

    time_model_start = time()

    bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])

    bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])

    for (i,bus) in ref[:bus]
        if length(ref[:bus_loads][i]) > 0
            bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
            bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
        end

        if length(ref[:bus_shunts][i]) > 0
            bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
            bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
        end
    end

    br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_ttm = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])

    for (i,branch) in ref[:branch]
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)

        br_g[i] = g
        br_b[i] = b

        br_tr[i] = tr
        br_ti[i] = ti
        br_ttm[i] = tr^2 + ti^2

        br_g_fr[i] = branch["g_fr"]
        br_b_fr[i] = branch["b_fr"]
        br_g_to[i] = branch["g_to"]
        br_b_to[i] = branch["b_to"]
    end

    var_lookup = Dict{String,Int}()

    var_init = Float64[]

    var_idx = 1
    for (i,bus) in ref[:bus]
        push!(var_init, 0.0) #va
        var_lookup["va_$(i)"] = var_idx
        var_idx += 1
        push!(var_init, 1.0) #vm
        var_lookup["vm_$(i)"] = var_idx
        var_idx += 1
    end

    for (i,gen) in ref[:gen]
        push!(var_init, 0.0) #pg
        var_lookup["pg_$(i)"] = var_idx
        var_idx += 1
        push!(var_init, 0.0) #qg
        var_lookup["qg_$(i)"] = var_idx
        var_idx += 1
    end

    for (l,i,j) in ref[:arcs]
        branch = ref[:branch][l]
        push!(var_init, 0.0) #p
        var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
        var_idx += 1
        push!(var_init, 0.0) #q
        var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
        var_idx += 1
    end

    @assert var_idx == length(var_init)+1

    function opf_objective(M, x)
        cost = 0.0
        for (i, gen) in ref[:gen]
            pg = x[var_lookup["pg_$(i)"]]
            cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
        end
        return cost
    end

    function opf_constraints_h(M, x)
        va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus])
        vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus])
        pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen])
        qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen])
        p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
        q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
        vm_fr = Dict(l => vm[branch["f_bus"]] for (l,branch) in ref[:branch])
        vm_to = Dict(l => vm[branch["t_bus"]] for (l,branch) in ref[:branch])
        va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch])
        va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch])
        return vcat(
            # @constraint(model, va[i] == 0)
            [va[i] for (i,bus) in ref[:ref_buses]],
            # @constraint(model,
            #     sum(p[a] for a in ref[:bus_arcs][i]) ==
            #     sum(pg[g] for g in ref[:bus_gens][i]) -
            #     sum(load["pd"] for load in bus_loads) -
            #     sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
            # )
            [
                sum(pg[j] for j in ref[:bus_gens][i]; init=0.0) -
                bus_pd[i] -
                bus_gs[i]*vm[i]^2 -
                sum(p[a] for a in ref[:bus_arcs][i])
                for (i,bus) in ref[:bus]
            ],
            # @constraint(model,
            #     sum(q[a] for a in ref[:bus_arcs][i]) ==
            #     sum(qg[g] for g in ref[:bus_gens][i]) -
            #     sum(load["qd"] for load in bus_loads) +
            #     sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
            # )
            [
                sum(qg[j] for j in ref[:bus_gens][i]; init=0.0) -
                bus_qd[i] +
                bus_bs[i]*vm[i]^2 -
                sum(q[a] for a in ref[:bus_arcs][i])
                for (i,bus) in ref[:bus]
            ],
            # @NLconstraint(model, p_fr ==  (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
            [
                (br_g[l]+br_g_fr[l])/br_ttm[l]*vm_fr[l]^2 +
                (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
                (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
                p[(l,i,j)]
                for (l,i,j) in ref[:arcs_from]
            ],
            # @NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
            [
                (br_g[l]+br_g_to[l])*vm_to[l]^2 +
                (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
                (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
                p[(l,i,j)]
                for (l,i,j) in ref[:arcs_to]
            ],
            # @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
            [
                -(br_b[l]+br_b_fr[l])/br_ttm[l]*vm_fr[l]^2 -
                (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
                (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
                q[(l,i,j)]
                for (l,i,j) in ref[:arcs_from]
            ],
            # @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
            [
                -(br_b[l]+br_b_to[l])*vm_to[l]^2 -
                (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
                (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
                q[(l,i,j)]
                for (l,i,j) in ref[:arcs_to]
            ],
        )
    end

    function opf_constraints_g(M, x)
        va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus])
        vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus])
        pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen])
        qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen])
        p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
        q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
        va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch])
        va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch])
        return vcat(
            [vm[i] - bus["vmax"] for (i,bus) in ref[:bus]],
            [bus["vmin"] - vm[i] for (i,bus) in ref[:bus]],
            [pg[i] - gen["pmax"] for (i,gen) in ref[:gen]],
            [gen["pmin"] - pg[i] for (i,gen) in ref[:gen]],
            [qg[i] - gen["qmax"] for (i,gen) in ref[:gen]],
            [gen["qmin"] - qg[i] for (i,gen) in ref[:gen]],
            [p[(l,i,j)] - ref[:branch][l]["rate_a"] for (l,i,j) in ref[:arcs]],
            [-ref[:branch][l]["rate_a"] - p[(l,i,j)] for (l,i,j) in ref[:arcs]],
            [q[(l,i,j)] - ref[:branch][l]["rate_a"] for (l,i,j) in ref[:arcs]],
            [-ref[:branch][l]["rate_a"] - q[(l,i,j)] for (l,i,j) in ref[:arcs]],
            # @constraint(model, va_fr - va_to <= branch["angmax"])
            [va_fr[l] - va_to[l] - ref[:branch][l]["angmax"] for (l,i,j) in ref[:arcs_from]],
            # @constraint(model, va_fr - va_to >= branch["angmin"])
            [ref[:branch][l]["angmin"] - (va_fr[l] - va_to[l]) for (l,i,j) in ref[:arcs_from]],
            # @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
            [p[(l,i,j)]^2 + q[(l,i,j)]^2 - ref[:branch][l]["rate_a"]^2 for (l,i,j) in ref[:arcs_from]],
            # @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
            [p[(l,i,j)]^2 + q[(l,i,j)]^2 - ref[:branch][l]["rate_a"]^2 for (l,i,j) in ref[:arcs_to]],
        )
    end

    model_variables = length(var_init)
    M = Euclidean(model_variables)

    model_constraints = length(opf_constraints_g(M, var_init)) +
                        length(opf_constraints_h(M, var_init))
    println("variables: $(model_variables)")
    println("constraints: $(model_constraints)")

    model_build_time = time() - time_model_start
    time_solve_start = time()

    solution = augmented_Lagrangian_method(
        M,
        opf_objective,
        (M, x) -> ForwardDiff.gradient(x -> opf_objective(M, x), x);
        g = opf_constraints_g,
        grad_g =  (M, x) -> ForwardDiff.jacobian(x -> opf_constraints_g(M, x), x),
        h = opf_constraints_h,
        grad_h =  (M, x) -> ForwardDiff.jacobian(x -> opf_constraints_h(M, x), x),
    )

    solve_time = time() - time_solve_start
    total_time = time() - time_data_start
    cost = opf_objective(M, solution)

    println("")
    println("\033[1mSummary\033[0m")
    println("   case........: $(file_name)")
    println("   variables...: $(model_variables)")
    println("   constraints.: $(model_constraints)")
    println("   feasible....: $(feasible)")
    println("   cost........: $(round(Int, cost))")
    println("   total time..: $(total_time)")
    println("     data time.: $(data_load_time)")
    println("     build time: $(model_build_time)")
    println("     solve time: $(solve_time)")
    # println("      callbacks: $(total_callback_time)")
    println("")

    return Dict(
        "case" => file_name,
        "variables" => model_variables,
        "constraints" => model_constraints,
        "feasible" => feasible,
        "cost" => cost,
        "time_total" => total_time,
        "time_data" => data_load_time,
        "time_build" => model_build_time,
        "time_solve" => solve_time,
        #"time_callbacks" => TBD,
    )
end

if isinteractive() == false
    solve_opf("$(@__DIR__)/data/pglib_opf_case5_pjm.m")
end

solve_opf("data/pglib_opf_case5_pjm.m")

[kellertuer]

I already wrote on Discourse such a code would have taken me weeks.

The main focus of Manopt is Optimisation on Manifolds, and sure constraint are possible but very rare (so rare that the first paper mentioning algorithms is to best of my knowledge from 2020). So the goal is different of that package.

AD is possible on manifolds but even there are more questions open. So sure one can do better if one is on Euclidean(n)

Constraints are – besides being new – also only in very few cases used. So that is also why all the fancy tools with upper and lower bounds are not available.

So then maybe it is not worth continuing, since the goal is too different.

But thanks for giving it a try on my very short start code :)

[ccoffrin]

This is a good exploration thank you @kellertuer and @odow!

@kellertuer, if you feel this application is out of scope for Manopt, we can close this topic. If you ever think it would be interesting to explore further we can revisit. I hope that rosetta-opf will have lasting value not only as a task for benchmarking but also as a way of understanding how different Julia optimization package work with respect to a common set of features (i.e., problem definition, derivative computation and solution reporting).

[kellertuer]

Thanks for your kind words.

I think the main problem is, that for now constraints are covered but not very performant – neither in setting the system up nor in how nice that is to set up. Due to the manifold structure I can – for now – also not say whether or when we will have anything beyond ALM and EPM.

Since I am also not knowledgable enough to fix why the benchmark does not finish (I even struggle to understand your constraints without any comments as the first post indicates), Manopt is maybe really not fitting so well here. It seems we just have a very small overlap (feature wise) where we just meet, but not well.

I have ManoptExampes.jl to hopefully illustrate how to use Manopt.jl in practical examples – besides that when we have the chance to get it running here that would be nice to contribute something to rosetta-ops for sure. We could go for just a few benchmarks maybe, since it does seem to work to some extend?

If you feel that is not enough and since I do not expect for Manopt to have more features in the constraint area in the next few years – unless a super-interesting topic with constraints or a very interested contributor comes along - we should maybe just close this for now.

edit: or to extend my response – I think we are a lot comparing apples to oranges here. If you want to solve the problems here and would ask me: “should we use Manopt.jl here?” My only answer would be: No. Nearly all other (Eculidean) solvers are better suited. They usually have a larger team of devs and much more methods fitting the problem here. Only if you get to the point where parts of your variables are a symmetric positive definite matrix or so, then consider Manopt.jl.

Otherwise it's a bit like taking an off-road electric bike and trying to go on an highway. Sure it can drive there. Will it be faster than the trucks and sports cars? Probably not. Do you still have to pedal (here AD maybe not as smooth as you think)? Sure, it‘s a bike. Is it uncomfortable with high speed (here many constraints)? Yes. But was it build for the highway? No. Should you go with an electric bike on an highway? Also probably no.

[kellertuer]

Looking at the code a bit, I think, one reason Manopt does not yet work well, is, that currently the gradients of the constraints are expected to be either

• a vector of gradients, one gradient per constraint
• a vector of functions, each of them providing a gradient for the corresponding constraint

The reason for that is again the manifold-setting – and not that I do not know what a Jacobian is. Gradients on manifolds might be represented in complicated ways, e.g. https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/fixedrankmatrices.html#Manifolds.UMVTVector, and it is not so easy to write that as a matrix actually.

So I think the correct form in the code above might be to split the Jacobian into a vector of its row vectors.

But maybe it is better to not continue this. In theory it should of course be possible, though Augmented Lagrangian is not the best choice (for the Euclidean world) for your problems. In practice the ALM in Manopt.jl has other things to fulfil which make it complicated to be adapted to the setting here.

[odow]

Closing for now. We can come back to this once https://github.com/JuliaManifolds/Manopt.jl/pull/264 is complete.