Removing P as a decision variable in the OPF problem

[gaadis]

I'm solving the OPF problem using the function _PMD.solve_mc_opf(networkfile, model, ipopt_solver). If I want to fix P i.e. remove it as a decision variable, so my optimization problem solves only Q, how do I go about it? Also, is there an objective function that minimizes losses or penalty function? The objective function that I see only minimizes cost. Any help in the right direction will be greatly appreciated.

[pseudocubic]

It sounds like you will want to write a custom problem, with some additional constraints and a new objective function. In order to not completely rewrite all the the existing constraints that call for active power variables, probably the best way would be to add some equality constraints to them. Here's an example of a modified version of the OPF problem (mostly cut and paste, with minor additions):

import PowerModelsDistribution as PMD
import JuMP
import Ipopt


function build_custom(pm::PMD.AbstractUnbalancedPowerModel)
    PMD.variable_mc_bus_voltage(pm)
    PMD.variable_mc_branch_power(pm)
    PMD.variable_mc_transformer_power(pm)
    PMD.variable_mc_switch_power(pm)
    PMD.variable_mc_generator_power(pm)
    PMD.variable_mc_load_power(pm)
    PMD.variable_mc_storage_power(pm)

    PMD.constraint_mc_model_voltage(pm)

    for i in PMD.ids(pm, :ref_buses)
        PMD.constraint_mc_theta_ref(pm, i)
    end

    # generators should be constrained before KCL, or Pd/Qd undefined
    for id in PMD.ids(pm, :gen)
        PMD.constraint_mc_generator_power(pm, id)
        constraint_generator_active_power_dispatch(pm, id)
    end

    # loads should be constrained before KCL, or Pd/Qd undefined
    for id in PMD.ids(pm, :load)
        PMD.constraint_mc_load_power(pm, id)
    end

    for i in PMD.ids(pm, :bus)
        PMD.constraint_mc_power_balance(pm, i)
    end

    for i in PMD.ids(pm, :storage)
        PMD.constraint_storage_state(pm, i)
        PMD.constraint_storage_complementarity_nl(pm, i)
        PMD.constraint_mc_storage_losses(pm, i)
        PMD.constraint_mc_storage_thermal_limit(pm, i)
        constraint_storage_active_power_dispatch(pm, id)
    end

    for i in PMD.ids(pm, :branch)
        PMD.constraint_mc_ohms_yt_from(pm, i)
        PMD.constraint_mc_ohms_yt_to(pm, i)

        PMD.constraint_mc_voltage_angle_difference(pm, i)

        PMD.constraint_mc_thermal_limit_from(pm, i)
        PMD.constraint_mc_thermal_limit_to(pm, i)
        PMD.constraint_mc_ampacity_from(pm, i)
        PMD.constraint_mc_ampacity_to(pm, i)
    end

    for i in PMD.ids(pm, :switch)
        PMD.constraint_mc_switch_state(pm, i)
        PMD.constraint_mc_switch_thermal_limit(pm, i)
        PMD.constraint_mc_switch_ampacity(pm, i)
    end

    for i in PMD.ids(pm, :transformer)
        PMD.constraint_mc_transformer_power(pm, i)
    end

    objective_custom_min(pm)
end

From this example, you would need three more functions, two to add the power constraints, and one for your objective.

function constraint_generator_active_power_dispatch(pm, i; nw=PMD.nw_id_default)
    for (idx,c) in enumerate(PMD.ref(pm, nw, :gen, i, "connections"))
        JuMP.@constraint(pm.model, PMD.var(pm, nw, :pg, i)[c] == PMD.ref(pm, nw, :gen, i, "pg")[idx])
    end
end


function constraint_storage_active_power_dispatch(pm, i; nw=PMD.nw_id_default)
    for (idx,c) in enumerate(PMD.ref(pm, nw, :storage, i, "connections"))
        JuMP.@constraint(pm.model, PMD.var(pm, nw, :ps, i)[c] == PMD.ref(pm, nw, :storage, i, "ps")[idx])
    end
end


function objective_custom_min(pm)
    JuMP.@objective(pm.model, Min, sum(sum(sum(qg) for (i, qg) in PMD.var(pm, n, :qg)) + sum(sum(qs) for (i, qs) in PMD.var(pm, n, :qs)) for n in PMD.nw_ids(pm)))
end

Here I have assumed that your desired active power setpoints are in the property "pg" on generators and "ps" on storage.

Your objective function is similarly straightforward. In this example I minimized the reactive powers from storage and generators.

Here is an example of this working, where I first find the power set points using OPF

pmd_path = joinpath(dirname(pathof(PMD)), "..")
data = PMD.parse_file("$pmd_path/test/data/opendss/case3_unbalanced.dss")

r_opf = PMD.solve_mc_opf(data, PMD.ACPUPowerModel, Ipopt.Optimizer)

for t in ["voltage_source", "generator", "solar"]
    for (i,obj) in get(r_opf["solution"], t, Dict())
        data[t][i]["pg"] = obj["pg"]
    end
end
for (i,obj) in get(r_opf["solution"], "storage", Dict())
    data["storage"][i]["ps"] = obj["ps"]
end

...then run your custom problem.

r_custom = PMD.solve_mc_model(data, PMD.ACPUPowerModel, Ipopt.Optimizer, build_custom)

One final note is that I found a minor bug in the ENGINEERING to MATHEMATICAL model transformation where voltage sources were not pulling in the pg/qg values. I will put a fix in today for that which will be required for my above example to work as-is.

[gaadis]

When I ran the code above, I got the error below:
ERROR: LoadError: UndefVarError: constraint_mc_ampacity_from not defined

I commented out the two constraints below,

    #PMD.constraint_mc_ampacity_from(pm, i)
    #PMD.constraint_mc_ampacity_to(pm, i)

and the code ran but converged to a point of local infeasibility as seen below

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.

I then, commented out the "constraint_generator_active_power_dispatch" function and the code converged with an optimal solution as seen below

Are you able to replicate the same issue I described above? How can I resolve this issue described above?

In a separate case, I used the objective function for the quadratic penalty for bus power slack variables as shown below;

but I got an undefined variable error with nws. Any idea on how to get rid of this error will be appreciated. Thank you.

[pseudocubic]

First, can you confirm which version of PMD you are importing?

The first error regarding constraint_mc_ampacity_from not found is definitely strange. I cannot replicate this.

The second error may just be because the pg / ps set-points on generation / storage are wrong?

The final error is likely just because if you are importing PMD like import PowerModelsDistribution as PMD then the PMD functions (var, ref, nws) should have PMD. in front of them.

[gaadis]

I'm using PMD v 0.11.4

[pseudocubic]

Is there a reason that you need to use that version, or is it possible to update to 0.14.9 (the latest release)?

Note that if you are stuck on that version, then the example function I wrote is not compatible and you may need to manually incorporate fixes that you are missing.

[gaadis]

When I install the PowerModelsDistributionStateEstimation package, I believe it downgrades PMD to v0.11.4. I'll update PMD to 0.14.9 and uninstall the PowerModelsDistributionStateEstimation to see what happens.

[gaadis]

I upgraded my PMD to 0.14.9 and it fixed the first problem about the ampacity. I'll change the setpoints to see if I can obtain feasible points. Also, when I use the objective function for the quadratic penalty for bus power slack variables as shown below.

I'm getting the following error even though I have PMD. in front of the function

[pseudocubic]

In your problem you will want to ensure that the function variable_mc_slack_bus_power gets called, e.g., https://github.com/lanl-ansi/PowerModelsDistribution.jl/blob/main/src/prob/debug.jl#L25

[gaadis]

That resolved the issue. Thanks.