Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Getting a feasible result for OPF function but an infeasible result with PF function using the same data #832

Closed
gaadis opened this issue Jul 29, 2022 · 7 comments

Comments

@gaadis
Copy link

gaadis commented Jul 29, 2022

I am running a 34-bus distribution dataset through the opf function, and I am getting a feasible result. However, when I run the same dataset through the pf function, the result is infeasible. I have combed through the data to figure out the problem but I'm not having any luck. Suggestions will be greatly appreciated. My code and data are posted below

###The error that I'm getting with the pf function #####
image

The feasible opf result with the same dataset

image

code

using JuMP, Ipopt
import PowerModels
import PowerModelsDistribution
const _PMs = PowerModels
const _PMD = PowerModelsDistribution

ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer,"max_cpu_time"=>300.0,
"tol"=>1e-10,
"print_level"=>0)
model = _PMD.ACPUPowerModel
pm_path = joinpath(dirname(pathof(PowerModels)), "..")
case5 = PowerModels.parse_file("$(pm_path)/test/data/matpower/case34Period1.m");
_PMD.make_multiconductor!(case5, 1)
pf_result = _PMD.solve_mc_pf(case5, model, ipopt_solver)
opf_result = _PMD.solve_mc_opf(case5, model, ipopt_solver)

Data set (case34Period1.m)

function mpc = case34Period1
mpc.version = '2';
mpc.baseMVA = 1;

%bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
mpc.gen = [
4 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
6 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
11 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
12 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
23 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
24 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
32 0.0000 0 10 -10 1 1 1 15.0 -5.6 0 0 0 0 0 0 0 0 0 0 0
];

%bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
mpc.bus = [
1 1 0.1537 0.0797 0 0 1 1 0.0 13.2 1 1.05 0.95
2 1 0.0926 0.0488 0 0 1 1 0.0 13.2 1 1.05 0.95
3 1 0.1020 0.0538 0 0 1 1 0.0 13.2 1 1.05 0.95
4 3 0.0313 0.0157 0 0 1 1 0.0 13.2 1 1.05 0.95
5 1 0.0268 0.0134 0 0 1 1 0.0 13.2 1 1.05 0.95
6 2 0.0000 0.0000 0 0 1 1 0.0 13.2 1 1.05 0.95
7 1 0.0000 0.0000 0 0 1 1 0.0 13.2 1 1.05 0.95
8 1 0.0000 0.0000 0 0 1 1 0.0 13.2 1 1.05 0.95
9 1 0.0095 0.0038 0 0 1 1 0.0 13.2 1 1.05 0.95
10 1 0.0457 0.0228 0 0 1 1 0.0 13.2 1 1.05 0.95
11 2 0.2763 0.1422 0 0 1 1 0.0 13.2 1 1.05 0.95
12 2 0.0000 0.0000 0 0 1 1 0.0 13.2 1 1.05 0.95
13 1 0.0865 0.0424 0 0 1 1 0.0 13.2 1 1.05 0.95
14 1 0.0578 0.0289 0 0 1 1 0.0 13.2 1 1.05 0.95
15 1 0.0154 0.0070 0 0 1 1 0.0 13.2 1 1.05 0.95
16 1 0.1443 0.0640 0 0 1 1 0.0 13.2 1 1.05 0.95
17 1 0.0082 0.0041 0 0 1 1 0.0 13.2 1 1.05 0.95
18 1 0.0061 0.0030 0 0 1 1 0.0 13.2 1 1.05 0.95
19 1 0.0000 0.0000 0 0 1 1 0.0 13.2 1 1.05 0.95
20 1 0.0259 0.0121 0 0 1 1 0.0 13.2 1 1.05 0.95
21 1 0.0000 0.0000 0 0 1 1 0.0 13.2 1 1.05 0.95
22 1 1.0973 0.5486 0 0 1 1 0.0 13.2 1 1.05 0.95
23 2 0.0755 0.0385 0 0 1 1 0.0 13.2 1 1.05 0.95
24 2 0.0036 0.0018 0 0 1 1 0.0 13.2 1 1.05 0.95
25 1 0.2732 0.1381 0 0 1 1 0.0 13.2 1 1.05 0.95
26 1 0.0128 0.0071 0 0 1 1 0.0 13.2 1 1.05 0.95
27 1 1.4573 1.1099 0 0 1 1 0.0 13.2 1 1.05 0.95
28 1 0.1326 0.0663 0 0 1 1 0.0 13.2 1 1.05 0.95
29 1 0.2144 0.1604 0 0 1 1 0.0 13.2 1 1.05 0.95
30 1 0.6136 0.3738 0 0 1 1 0.0 13.2 1 1.05 0.95
31 1 0.2216 0.1144 0 0 1 1 0.0 13.2 1 1.05 0.95
32 2 0.1785 0.1177 0 0 1 1 0.0 13.2 1 1.05 0.95
33 1 0.0473 0.0236 0 0 1 1 0.0 13.2 1 1.05 0.95
34 1 0.0500 0.0250 0 0 1 1 0.0 13.2 1 1.05 0.95
];

%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [ %% (r, x and b specified in p.u. here)
1 2 0.000107 0.000102 0.015709 0 0 0 0 0 1 -60 60;
2 3 0.000072 0.000068 0.010533 0 0 0 0 0 1 -60 60;
3 4 0.001340 0.001273 0.196240 0 0 0 0 0 1 -60 60;
4 5 0.000496 0.000263 0.028796 0 0 0 0 0 1 -60 60;
4 6 0.001559 0.001481 0.228320 0 0 0 0 0 1 -60 60;
6 7 0.001236 0.001174 0.181010 0 0 0 0 0 1 -60 60;
9 10 0.000146 0.000078 0.0084839 0 0 0 0 0 1 -60 60;
9 13 0.000610 0.000431 0.059795 0 0 0 0 0 1 -60 60;
10 11 0.004118 0.002185 0.238890 0 0 0 0 0 1 -60 60;
11 12 0.001175 0.000623 0.068169 0 0 0 0 0 1 -60 60;
13 14 0.000259 0.000137 0.015033 0 0 0 0 0 1 -60 60;
13 15 0.000050 0.000035 0.0049194 0 0 0 0 0 1 -60 60;
15 16 0.001222 0.000862 0.119701 0 0 0 0 0 1 -60 60;
16 17 0.000031 0.000022 0.0030454 0 0 0 0 0 1 -60 60;
20 23 0.000293 0.000207 0.028697 0 0 0 0 0 1 -60 60;
25 30 0.000121 0.000085 0.011830 0 0 0 0 0 1 -60 60;
25 26 0.000017 0.000012 0.0016398 0 0 0 0 0 1 -60 60;
31 32 0.000051 0.000036 0.0050366 0 0 0 0 0 1 -60 60;
31 33 0.000017 0.000012 0.0016398 0 0 0 0 0 1 -60 60;
26 27 0.000081 0.000057 0.0079062 0 0 0 0 0 1 -60 60;
27 28 0.000218 0.000154 0.021318 0 0 0 0 0 1 -60 60;
28 29 0.000032 0.000022 0.0031039 0 0 0 0 0 1 -60 60;
8 9 0.000019 0.000013 0.0018155 0 0 0 0 0 1 -60 60;
17 18 0.001995 0.001059 0.115750 0 0 0 0 0 1 -60 60;
17 19 0.002201 0.001554 0.215690 0 0 0 0 0 1 -60 60;
23 24 0.000139 0.000074 0.0080374 0 0 0 0 0 1 -60 60;
23 25 0.000348 0.000246 0.034143 0 0 0 0 0 1 -60 60;
30 31 0.000160 0.000113 0.015695 0 0 0 0 0 1 -60 60;
33 34 0.000285 0.000211 0.024954 0 0 0 0 0 1 -60 60;
21 22 0.015729 0.014939 0.0017946 0 0 0 0 0 1 -60 60;
20 21 0.003800 0.008160 0 0 0 0 0 0 1 -60 60;
];

%% generator cost data
% 1 startup shutdown n x1 y1 ... xn yn
% 2 startup shutdown n c(n-1) ... c0
mpc.gencost = [
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
2 0.0 0.0 4 0.000100000 0.00100000 0.0100000 0.100000;
% if n =4, then (0.1 x 10^3) + (4 x 10^2) + (0.0 x 10^1) + (2 x 10^0)
];

@ccoffrin
Copy link
Member

Hi @gaadis, should this issue should be moved to PowerModelsDistribution? It looks to me that is the package you are using in this example.

At a high level, what you observe is reasonable and common in many datasets. The OPF problem provides more degrees of freedom to find a feasible power flow. So if the generator injection set points or voltage profile of your dataset are not setup correctly, AC-PF will fail but AC-OPF will work.

CC @pseudocubic

@gaadis
Copy link
Author

gaadis commented Jul 29, 2022

Yes, please move it to PMD

@gaadis
Copy link
Author

gaadis commented Jul 29, 2022

Hi @gaadis, should this issue should be moved to PowerModelsDistribution? It looks to me that is the package you are using in this example.

At a high level, what you observe is reasonable and common in many datasets. The OPF problem provides more degrees of freedom to find a feasible power flow. So if the generator injection set points or voltage profile of your dataset are not setup correctly, AC-PF will fail but AC-OPF will work.

CC @pseudocubic

What will be the best way to check if the generator injector set points or voltage profile of the dataset is setup properly?

@ccoffrin
Copy link
Member

Usually the best way to check is exactly how you have done it here. Run an AC-PF, if it fails, try an AC-OPF, if that works then most likely the set-point in your AC-PF data is not suitable.

@gaadis
Copy link
Author

gaadis commented Jul 29, 2022

Usually the best way to check is exactly how you have done it here. Run an AC-PF, if it fails, try an AC-OPF, if that works then most likely the set-point in your AC-PF data is not suitable.

Is there an easier way to determine the right set-points instead of trial and error?

@pseudocubic
Copy link
Collaborator

pseudocubic commented Aug 1, 2022

It looks like you are using a workflow that isn't intended to be well supported in PowerModelsDistribution. If you note the documentation for make_multiconductor!, it says that it is a "Hacky helper function" to convert into multiconductor data. Most notably, it of course doesn't include any off-diagonal terms in the impedance. I will make a note to make the unsupported nature of this function more clear in the documentation.

However, from my understanding of the code you posted, you are converting a matpower input file, parsed with PowerModels, into a multiconductor format with only one conductor. Therefore, running PMD.solve_mc_pf should be equivalent to PM.solve_pf in this case anyway, where ACPUPowerModel PowerModelsDistribution in is equivalent to ACPPowerModel in PowerModels.

When I attempt to run the network model in PowerModels (without converting to multiconductor, i.e., case5 = PowerModels.parse_file("case34period1.m");), I similarly do not get a result from solve_pf:

julia> pm_pf_result = _PMs.solve_pf(case5, _PMs.ACPPowerModel, ipopt_solver)
Dict{String, Any} with 8 entries:
  "solve_time"         => 10.0657
  "optimizer"          => "Ipopt"
  "termination_status" => ITERATION_LIMIT
  "dual_status"        => UNKNOWN_RESULT_STATUS
  "primal_status"      => UNKNOWN_RESULT_STATUS
  "objective"          => 0.0
  "solution"           => Dict{String, Any}("baseMVA"=>1, "gen"=>Dict{String, Any}("4"=>Dict{String, Any}("qg"=>-0.0316291, "pg"=>8.76596e-8), 
  "objective_lb"       => -Inf

but get this from AC-OPF:

julia> pm_opf_result = _PMs.solve_opf(case5, _PMs.ACPPowerModel, ipopt_solver)
Dict{String, Any} with 8 entries:
  "solve_time"         => 0.25114
  "optimizer"          => "Ipopt"
  "termination_status" => LOCALLY_SOLVED
  "dual_status"        => FEASIBLE_POINT
  "primal_status"      => FEASIBLE_POINT
  "objective"          => 0.763219
  "solution"           => Dict{String, Any}("baseMVA"=>1, "branch"=>Dict{String, Any}("24"=>Dict{String, Any}("qf"=>-0.101484, "qt"=>-0.003, "p…
  "objective_lb"       => -Inf

When I try to update the setpoints from the OPF result though, I get a numerical error:

julia> for (i,gen) in pm_opf_result["solution"]["gen"]
           case5["gen"][i]["pg"] = gen["pg"]
           case5["gen"][i]["qg"] = gen["qg"]
       end

julia> pm_pf_result = _PMs.solve_pf(case5, _PMs.ACPPowerModel, ipopt_solver)
Dict{String, Any} with 8 entries:
  "solve_time"         => 2.39664
  "optimizer"          => "Ipopt"
  "termination_status" => NUMERICAL_ERROR
  "dual_status"        => UNKNOWN_RESULT_STATUS
  "primal_status"      => UNKNOWN_RESULT_STATUS
  "objective"          => 0.0
  "solution"           => Dict{String, Any}("baseMVA"=>1, "gen"=>Dict{String, Any}("4"=>Dict{String, Any}("qg"=>-0.644033, "pg"=>0.323762), "1"
  "objective_lb"       => -Inf

pseudocubic added a commit to lanl-ansi/PowerModelsDistribution.jl that referenced this issue Aug 1, 2022
Makes docstring of make_multinetwork! more clear that this is an unsupported avenue to create valid multiconductor networks

Relates to [PowerModels.jl#832](lanl-ansi/PowerModels.jl#832)
@ccoffrin
Copy link
Member

Closing due to no actionable items for PowerModels other than #872.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants