diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 62b09692..56d88b00 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,7 +13,7 @@ jobs: fail-fast: false matrix: version: - - '1.8' + - '1.9' - '1.6' - '1' os: diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index d2617481..cad05b03 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -1,26 +1,26 @@ name: Documentation + on: push: - branches: [master] + branches: + - master tags: '*' pull_request: - types: [opened, synchronize, reopened] - schedule: - - cron: '0,5,10 13 * * *' + jobs: - docs: - if: "!contains(github.event.head_commit.message, 'skip ci')" - name: Documentation + build: + permissions: + contents: write + statuses: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@latest + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 with: version: '1' - - name: Install Dependencies - run: julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - - name: Build and Deploy + - name: Install dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} - run: julia --project=docs docs/make.jl \ No newline at end of file + run: julia --project=docs/ docs/make.jl diff --git a/Project.toml b/Project.toml index 6c79650c..e5423af2 100644 --- a/Project.toml +++ b/Project.toml @@ -2,13 +2,15 @@ name = "PowerModelsACDC" uuid = "ff45984e-d068-5f4c-9e32-c4133509d236" autors = ["Hakan Ergun", "Frederik Geth", "Jay Dave"] repo = "https://github.com/Electa-Git/PowerModelsACDC.jl" -version = "0.6.3" +version = "0.7.0" [deps] InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Memento = "f28f55f0-a522-5efc-85c2-fe41dfb9b2d9" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] HiGHS = "~0.3, ~1.4" @@ -17,8 +19,10 @@ Ipopt = "~0.8, ~0.9, 1" JuMP = "~0.22.0, ~0.23, ~1" Juniper = "~0.8, ~0.9" Memento = "~1.0, ~1.1, ~1.2, ~1.3, ~1.4" +NLsolve = "~4.5" PowerModels = "~0.19" julia = "^1" +SparseArrays = "1" [extras] Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" @@ -28,4 +32,4 @@ Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Ipopt", "HiGHS", "Juniper"] \ No newline at end of file +test = ["Test", "Ipopt", "HiGHS", "Juniper"] diff --git a/README.md b/README.md index c6af1b99..5d2f2b19 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ PowerModelsACDC.jl is a Julia/JuMP/PowerModels package with models for DC lines, Building upon the PowerModels architecture, the code is engineered to decouple problem specifications (e.g. Power Flow, Optimal Power Flow, ...) from the power network formulations (e.g. AC, DC-approximation, SOC-relaxation, ...). **Installation** -The latest stable release of PowerModelACDC can be installed using the Julia package manager with +The latest stable release of PowerModelsACDC can be installed using the Julia package manager with ```julia Pkg.add("PowerModelsACDC") @@ -56,7 +56,7 @@ The developers thank Carleton Coffrin (LANL) for his support. - Hakan Ergun (KU Leuven / EnergyVille): Main developer - Frederik Geth (KU Leuven / EnergyVille): Formulations & relaxations of the OPF problem - Jay Dave (KU Leuven / EnergyVille): Transmission expansion plannning -- Ghulam Mohy Ud Din (CSIRO): ACR formulation of the OPF problem +- Ghulam Mohy Ud Din (CSIRO): ACR formulation of the OPF problem, Sequential AC/DC Grid Power Flow using NLsolve ## Citing PowerModelsACDC @@ -89,7 +89,7 @@ month={July},} ISSN = {1751-8687}, title = {TNEP of meshed HVDC grids: ‘AC’, ‘DC’ and convex formulations}, journal = {IET Generation, Transmission & Distribution}, - issue = {24}, + issue = {24}, volume = {13}, year = {2019}, month = {December}, diff --git a/docs/.gitignore b/docs/.gitignore deleted file mode 100644 index e69de29b..00000000 diff --git a/docs/README.md b/docs/README.md index db8cafbf..31f67d1d 100644 --- a/docs/README.md +++ b/docs/README.md @@ -1,21 +1,36 @@ -# Building the Documentation for PowerModels.jl +# Documentation for PowerModelsACDC.jl -## Installation -We rely on [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl). To install it, run the following command in a julia session: +You can read this documentation online at +. -```julia -Pkg.add("Documenter") -``` +## Preview the documentation (for developers) -## Building the Docs -To preview the html output of the documents, run the following command: +While developing PowerModelsACDC.jl, you can preview the documentation locally in your +browser with live-reload capability. +In other words, when you modify a file, every browser tab that is currently displaying the +corresponding page is automatically refreshed. -```julia -julia --color=yes make.jl -``` +### Instructions for *nix -You can then view the documents in `build/index.html`. +1. Copy the following zsh/Julia code snippet: -**Warning**: Do not `git commit` the contents of build (or any other content generated by Documenter) to your repository's master branch. This helps to avoid including unnessesary changes for anyone reviewing commits that happen to include documentation changes. + ```julia + #!/bin/zsh + #= # The following line is zsh code + julia -i $0:a # The string `$0:a` represents this file in zsh + =# # Following lines are Julia code + import Pkg + Pkg.activate(; temp=true) + Pkg.develop("PowerModelsACDC") + Pkg.add("Documenter") + Pkg.add("LiveServer") + using PowerModelsACDC, LiveServer + cd(dirname(dirname(pathof(PowerModelsACDC)))) + servedocs() + exit() + ``` -For further details, please read the [documentation for Documenter.jl](https://juliadocs.github.io/Documenter.jl/stable/). +2. Save it as a zsh script (name it like `preview_powermodelsacdc_docs.sh`). +3. Assign execute permission to the script: `chmod u+x preview_powermodelsacdc_docs.sh`. +4. Run the script. +5. Open your favorite web browser and navigate to `http://localhost:8000`. diff --git a/docs/_config.yml b/docs/_config.yml deleted file mode 100644 index c7418817..00000000 --- a/docs/_config.yml +++ /dev/null @@ -1 +0,0 @@ -theme: jekyll-theme-slate \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index a8d4363b..222bd0c2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,30 +1,29 @@ using Documenter, PowerModelsACDC -Documenter.makedocs( - modules = PowerModelsACDC, - format = Documenter.HTML(), +makedocs( + modules = [PowerModelsACDC], sitename = "PowerModelsACDC", - authors = "Frederik Geth, Hakan Ergun", + warnonly = :missing_docs, pages = [ - "Home" => "index.md", + "Home" => "index.md" "Manual" => [ - "Getting Started" => "quickguide.md", - "Results" => "result-data.md", - ], + "Getting Started" => "quickguide.md" + "Results" => "result-data.md" + ] "Library" => [ - "Network Formulations" => "formulations.md", - "Problem Specifications" => "specifications.md", - "Problem Types" => "problems.md", + "Network Formulations" => "formulations.md" + "Problem Specifications" => "specifications.md" + "Problem Types" => "problems.md" "Modeling Components" => [ - "Objective" => "objective.md", - "Variables" => "variables.md", - "Constraints" => "constraints.md", - ], + "Objective" => "objective.md" + "Variables" => "variables.md" + "Constraints" => "constraints.md" + ] "File IO" => "parser.md" - ], + ] ] ) -Documenter.deploydocs( +deploydocs( repo = "github.com/Electa-Git/PowerModelsACDC.jl.git" ) diff --git a/docs/src/constraints.md b/docs/src/constraints.md index 7f47e72a..9620ddb6 100644 --- a/docs/src/constraints.md +++ b/docs/src/constraints.md @@ -5,12 +5,6 @@ All the OPF constraints for the AC grids have been re-used from PowerModels.jl, CurrentModule = PowerModelsACDC ``` -## Unit Constraints - -```@docs -constraint_active_load_gen_aggregation -``` - ## DC Bus Constraints ### Setpoint Constraints diff --git a/docs/src/formulations.md b/docs/src/formulations.md index ae507f4f..c77c7b7d 100644 --- a/docs/src/formulations.md +++ b/docs/src/formulations.md @@ -60,7 +60,7 @@ $Q^{conv, ac} = \sin\varphi_{c} \cdot S^{conv,ac,rated}$ ### ACDC converters Two separate current variables, $I^{conv, ac}$ and $i^{conv, ac, sq}$ are defined, the nonconvex relation $i^{conv, ac, sq} = (I^{conv, ac})^2$ is convexified. -- Linking both current variables: $(I^{conv, ac})^2$ $\leq$ $i^{conv, ac, sq}$ +- Linking both current variables: $(I^{conv, ac})^2$ $=$ $i^{conv, ac, sq}$ - Power balance: $P^{conv, ac}_{ij} + P^{conv, dc}_{ji}$ = $a + b\cdot I^{conv, ac} + c\cdot i^{conv, ac, sq}$. - Converter current variable model: $(P^{conv,ac}_{ij})^2$ + $(Q^{conv,ac}_{ij})^2$ = $(U_{ri}^2+U_{ii}^2) \cdot i^{conv, ac, sq}$. - LCC converters, active /reactive power: Same model as ACP formulation diff --git a/docs/src/index.md b/docs/src/index.md index 4ad01281..2afd3dfc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,4 +1,4 @@ -# PowerModelACDC.jl Documentation +# PowerModelsACDC.jl Documentation ```@meta CurrentModule = PowerModelsACDC @@ -26,16 +26,17 @@ Developed by: - Hakan Ergun, Jay Dave, KU Leuven / EnergyVille - Frederik Geth, CSIRO -## Installation of PowerModelACDC +## Installation of PowerModelsACDC -The latest stable release of PowerModelACDC can be installed using the Julia package manager with +The latest stable release of PowerModelsACDC can be installed using the Julia package manager with ```julia Pkg.add("PowerModelsACDC") ``` -!!! Important +!!! note This is a research-grade optimization package. ## Special Thanks To -Jef Beerten (KU Leuven/EnergyVille) for his insights in AC/DC power flow modelling. -Carleton Coffrin (Los Alamos National Laboratory) for his countless design tips. + +- Jef Beerten (KU Leuven/EnergyVille), for his insights in AC/DC power flow modelling. +- Carleton Coffrin (Los Alamos National Laboratory), for his countless design tips. diff --git a/docs/src/objective.md b/docs/src/objective.md index 2c64f403..cb530e6f 100644 --- a/docs/src/objective.md +++ b/docs/src/objective.md @@ -9,12 +9,18 @@ * Minimisation of consturction cost for HVDC branches and converters and generation costs and AC branches(acdc_tnep) - ```@meta CurrentModule = PowerModelsACDC ``` +```@docs +objective_min_fuel_cost +``` + +```@docs +objective_min_cost +``` ```@docs -_PM.objective_min_fuel_cost +objective_min_cost_acdc ``` diff --git a/docs/src/problems.md b/docs/src/problems.md index e304aad0..586293d9 100644 --- a/docs/src/problems.md +++ b/docs/src/problems.md @@ -11,6 +11,11 @@ run_acdcopf(file or data, formulation, solver, setting) run_acdcpf(file or data, formulation, solver, setting) ``` +## Sequential Hybrid AC/DC power flow (Native) +```julia +run_sacdcpf(file or data) +``` + ## DC grid TNEP problem (optimal placement of converters and dc branches) ```julia run_tnepopf(file or data, formulation, solver, setting) diff --git a/docs/src/specifications.md b/docs/src/specifications.md index cd4c4fab..f2c9686f 100644 --- a/docs/src/specifications.md +++ b/docs/src/specifications.md @@ -1,5 +1,35 @@ # Problem Specifications + +## ACDCPF + +PF with support for AC and DC grids at the same time, including AC/DC converters. + +## Generic AC DC Power Flow + +The general purpose ac dc power flow solver in PowerModelsACDC is, + +```@docs +run_acdcpf +``` + +This function builds a JuMP model for a wide variety of the power flow formulations +supported by PowerModelsACDC. + +## Sequential AC DC Power Flow (Native) + +The sequential ac dc power flow solver in PowerModelsACDC uses the package +[NLSolve](https://github.com/JuliaNLSolvers/NLsolve.jl) for solving the AC +DC power flow problem in `ACPPowerModel` formulation sequentially. + +```@docs +run_sacdcpf +``` + +!!! tip + If `run_sacdcpf` fails to converge try `run_acdcpf` instead. + + ## ACDCOPF OPF with support for AC and DC grids at the same time, including AC/DC converters. diff --git a/src/PowerModelsACDC.jl b/src/PowerModelsACDC.jl index a1314c77..de24e11d 100644 --- a/src/PowerModelsACDC.jl +++ b/src/PowerModelsACDC.jl @@ -10,9 +10,11 @@ const _PM = PowerModels import InfrastructureModels # import InfrastructureModels: ids, ref, var, con, sol, nw_ids, nws, optimize_model!, @im_fields const _IM = InfrastructureModels +import SparseArrays +import NLsolve -import JuMP: with_optimizer -export with_optimizer +import JuMP: optimizer_with_attributes +export optimizer_with_attributes # Create our module level logger (this will get precompiled) const _LOGGER = Memento.getlogger(@__MODULE__) @@ -30,6 +32,7 @@ include("prob/tnepopf.jl") include("prob/tnepopf_bf.jl") include("prob/mp_tnepopf.jl") include("prob/mp_tnepopf_bf.jl") +include("prob/sacdcpf.jl") diff --git a/src/core/base.jl b/src/core/base.jl index 877e073a..0861e573 100644 --- a/src/core/base.jl +++ b/src/core/base.jl @@ -134,7 +134,6 @@ function find_all_ac_grids(branches_ac, buses_ac) ACgrids["1"]["Buses"] = [branches_ac[1]["f_bus"] branches_ac[1]["t_bus"]] closed_buses = [branches_ac[1]["f_bus"] branches_ac[1]["t_bus"]] closed_branches = [1] - connections = [] buses = [] for (i, bus) in buses_ac if VERSION < v"0.7.0-" diff --git a/src/core/constraint.jl b/src/core/constraint.jl index 233df875..1da42c14 100644 --- a/src/core/constraint.jl +++ b/src/core/constraint.jl @@ -34,9 +34,19 @@ end ###################### TNEP Constraints ############################ +"do nothing, this model does not have complex voltage constraints" function constraint_voltage_dc_ne(pm::_PM.AbstractPowerModel, n::Int) end +""" +``` +z[c] * lb[c] <= pconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_pf_fr[c] <= z[c] * ub[c] +``` +""" function constraint_converter_limit_on_off(pm::_PM.AbstractDCPModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax) pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i] pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i] @@ -56,7 +66,21 @@ function constraint_converter_limit_on_off(pm::_PM.AbstractDCPModel, n::Int, i, JuMP.@constraint(pm.model, pconv_pr_fr >= pmin * z) end - +""" +``` +z[c] * lb[c] <= pconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= qconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_pf_fr[c] <= z[c] * ub[c] +``` +""" function constraint_converter_limit_on_off(pm::_PM.AbstractACPModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax) pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i] pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i] @@ -97,7 +121,24 @@ function constraint_converter_limit_on_off(pm::_PM.AbstractACPModel, n::Int, i, end - +""" +``` +z[c] * lb[c] <= pconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= qconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= iconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c]^2 <= iconv_sq[c] <= z[c] * ub[c]^2 +``` +""" function constraint_converter_limit_on_off(pm::_PM.AbstractBFModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax) #converter pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i] @@ -150,7 +191,24 @@ function constraint_converter_limit_on_off(pm::_PM.AbstractBFModel, n::Int, i, p JuMP.@constraint(pm.model, ipr <= (bigM * imax)^2 * z) #big M = 2 end - +""" +``` +z[c] * lb[c] <= pconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= qconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= iconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c]^2 <= iconv_sq[c] <= z[c] * ub[c]^2 +``` +""" function constraint_converter_limit_on_off(pm::_PM.AbstractWModels, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax) #converter pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i] @@ -191,7 +249,24 @@ function constraint_converter_limit_on_off(pm::_PM.AbstractWModels, n::Int, i, p JuMP.@constraint(pm.model, qconv_pr_fr >= qmin * z) end - +""" +``` +z[c] * lb[c] <= pconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= qconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= iconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c]^2 <= iconv_sq[c] <= z[c] * ub[c]^2 +``` +""" function constraint_converter_limit_on_off(pm::_PM.AbstractWRMModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax) #converter pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i] @@ -232,6 +307,23 @@ function constraint_converter_limit_on_off(pm::_PM.AbstractWRMModel, n::Int, i, JuMP.@constraint(pm.model, qconv_pr_fr >= qmin * z) end +""" +``` +z[c] * lb[c] <= pconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= pconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= qconv_ac[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_dc[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_fr[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_tf_to[c] <= z[c] * ub[c] +z[c] * lb[c] <= qconv_pf_fr[c] <= z[c] * ub[c] + +z[c] * lb[c] <= iconv_ac[c] <= z[c] * ub[c] +``` +""" function constraint_converter_limit_on_off(pm::_PM.AbstractLPACModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax) #converter pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i] @@ -270,7 +362,13 @@ function constraint_converter_limit_on_off(pm::_PM.AbstractLPACModel, n::Int, i, JuMP.@constraint(pm.model, qconv_pr_fr >= qmin * z) end - +""" +``` +z[b] * lb[b] <= pfr[b] <= z[b] * ub[b] +z[b] * lb[b] <= pto[b] <= z[b] * ub[b] +z[b] * lb[b] <= ccm[b] <= z[b] * ub[b] +``` +""" function constraint_branch_limit_on_off(pm::_PM.AbstractBFModel, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin) p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx] p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx] @@ -283,7 +381,13 @@ function constraint_branch_limit_on_off(pm::_PM.AbstractBFModel, n::Int, i, f_id JuMP.@constraint(pm.model, ccm_dcgrid <= imax^2 * z) JuMP.@constraint(pm.model, ccm_dcgrid >= imin^2 * z) - end +end +""" +``` +z[b] * lb[b] <= pfr[b] <= z[b] * ub[b] +z[b] * lb[b] <= pto[b] <= z[b] * ub[b] +``` +""" function constraint_branch_limit_on_off(pm::_PM.AbstractWRModels, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin) p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx] p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx] @@ -293,7 +397,12 @@ function constraint_branch_limit_on_off(pm::_PM.AbstractWRModels, n::Int, i, f_i JuMP.@constraint(pm.model, p_to <= pmax * z) JuMP.@constraint(pm.model, p_to >= pmin * z) end - +""" +``` +z[b] * lb[b] <= pfr[b] <= z[b] * ub[b] +z[b] * lb[b] <= pto[b] <= z[b] * ub[b] +``` +""" function constraint_branch_limit_on_off(pm::_PM.AbstractPowerModel, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin) p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx] p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx] @@ -304,7 +413,12 @@ function constraint_branch_limit_on_off(pm::_PM.AbstractPowerModel, n::Int, i, f JuMP.@constraint(pm.model, p_to <= pmax * z) JuMP.@constraint(pm.model, p_to >= pmin * z) end - +""" +``` +z[b] * lb[b] <= pfr[b] <= z[b] * ub[b] +z[b] * lb[b] <= pto[b] <= z[b] * ub[b] +``` +""" function constraint_branch_limit_on_off(pm::_PM.AbstractACPModel, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin) p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx] p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx] @@ -316,7 +430,11 @@ function constraint_branch_limit_on_off(pm::_PM.AbstractACPModel, n::Int, i, f_i JuMP.@constraint(pm.model, p_to >= pmin * z) end - +""" +``` +z[c,n] = z[c,n-1] ∀ n > 1 +``` +""" function constraint_candidate_converters_mp(pm::_PM.AbstractPowerModel, n::Int, i::Int) z = _PM.var(pm, n, :conv_ne, i) z_1 = _PM.var(pm, n-1, :conv_ne, i) @@ -324,6 +442,11 @@ function constraint_candidate_converters_mp(pm::_PM.AbstractPowerModel, n::Int, JuMP.@constraint(pm.model, z == z_1) end +""" +``` +z[b,n] = z[b,n-1] ∀ n > 1 +``` +""" function constraint_candidate_dcbranches_mp(pm::_PM.AbstractPowerModel, n::Int, i::Int) z = _PM.var(pm, n, :branchdc_ne, i) z_1 = _PM.var(pm, n-1, :branchdc_ne, i) @@ -331,6 +454,11 @@ function constraint_candidate_dcbranches_mp(pm::_PM.AbstractPowerModel, n::Int, JuMP.@constraint(pm.model, z == z_1) end +""" +``` +z[b,n] = z[b,n-1] ∀ n > 1 +``` +""" function constraint_candidate_acbranches_mp(pm::_PM.AbstractPowerModel, n::Int, i::Int) z = _PM.var(pm, n, :branch_ne, i) z_1 = _PM.var(pm, n-1, :branch_ne, i) @@ -338,6 +466,11 @@ function constraint_candidate_acbranches_mp(pm::_PM.AbstractPowerModel, n::Int, JuMP.@constraint(pm.model, z == z_1) end +""" +``` +sum(p_dcgrid[a] for a in bus_arcs_dcgrid) + sum(p_dcgrid_ne[a] for a in bus_arcs_dcgrid_ne) + sum(pconv_dc[c] for c in bus_convs_dc) + sum(pconv_dc_ne[c] for c in bus_convs_dc_ne) == (-pd) +``` +""" function constraint_power_balance_dc_dcne(pm::_PM.AbstractPowerModel, n::Int, i::Int, bus_arcs_dcgrid, bus_arcs_dcgrid_ne, bus_convs_dc, bus_convs_dc_ne, pd) p_dcgrid = _PM.var(pm, n, :p_dcgrid) p_dcgrid_ne = _PM.var(pm, n, :p_dcgrid_ne) @@ -351,7 +484,11 @@ function constraint_power_balance_dc_dcne(pm::_PM.AbstractPowerModel, n::Int, i: end end - +""" +``` +sum(p_dcgrid_ne[a] for a in bus_arcs_dcgrid_ne) + sum(pconv_dc_ne[c] for c in bus_ne_convs_dc_ne) == (-pd_ne) +``` +""" function constraint_power_balance_dcne_dcne(pm::_PM.AbstractPowerModel, n::Int, i::Int, bus_arcs_dcgrid_ne, bus_ne_convs_dc_ne, pd_ne) p_dcgrid_ne = _PM.var(pm, n, :p_dcgrid_ne) pconv_dc_ne = _PM.var(pm, n, :pconv_dc_ne) diff --git a/src/core/objective.jl b/src/core/objective.jl index 10fbb582..f6a4b6b0 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -1,4 +1,4 @@ -"" +" Minimum fuel cost objective using lienar and quadratic terms or piecewise linear functions, as defined in Matpower" function objective_min_fuel_cost(pm::_PM.AbstractPowerModel) model = _PM.check_cost_models(pm) if model == 1 @@ -11,7 +11,6 @@ function objective_min_fuel_cost(pm::_PM.AbstractPowerModel) end -"" function objective_min_polynomial_fuel_cost(pm::_PM.AbstractPowerModel) order = _PM.calc_max_cost_index(pm.data)-1 @@ -23,7 +22,6 @@ function objective_min_polynomial_fuel_cost(pm::_PM.AbstractPowerModel) error("cost model order of $(order) is not supported") end end - function _objective_min_polynomial_fuel_cost_linear(pm::_PM.AbstractPowerModel) from_idx = Dict() for (n, nw_ref) in nws(pm) @@ -37,7 +35,6 @@ function _objective_min_polynomial_fuel_cost_linear(pm::_PM.AbstractPowerModel) for (n, nw_ref) in nws(pm)) ) end -"" function _objective_min_polynomial_fuel_cost_quadratic(pm::_PM.AbstractPowerModel) from_idx = Dict() for (n, nw_ref) in nws(pm) @@ -52,7 +49,6 @@ function _objective_min_polynomial_fuel_cost_quadratic(pm::_PM.AbstractPowerMode for (n, nw_ref) in nws(pm)) ) end -"" function objective_min_polynomial_fuel_cost(pm::_PM.AbstractConicModel) _PM.check_polynomial_cost_models(pm) pg_sqr = Dict() @@ -77,7 +73,6 @@ function objective_min_polynomial_fuel_cost(pm::_PM.AbstractConicModel) for (n, nw_ref) in _PM.nws(pm)) ) end -"" function objective_min_pwl_fuel_cost(pm::_PM.AbstractPowerModel) for (n, nw_ref) in _PM.nws(pm) @@ -104,6 +99,7 @@ function objective_min_pwl_fuel_cost(pm::_PM.AbstractPowerModel) end ##################### TNEP Objective ################### +" Objective consisting of generation cost + investment costs of HVDC lines and converters" function objective_min_cost(pm::_PM.AbstractPowerModel) gen_cost = Dict() for (n, nw_ref) in _PM.nws(pm) @@ -134,6 +130,7 @@ function objective_min_cost(pm::_PM.AbstractPowerModel) ) end +" Objective consisting of generation cost + investment costs of HVDC lines and converters + investment costs of AC lines" function objective_min_cost_acdc(pm::_PM.AbstractPowerModel) gen_cost = Dict() for (n, nw_ref) in _PM.nws(pm) diff --git a/src/core/relaxation_scheme.jl b/src/core/relaxation_scheme.jl index 2ca0da62..6d4da176 100644 --- a/src/core/relaxation_scheme.jl +++ b/src/core/relaxation_scheme.jl @@ -9,7 +9,7 @@ function relaxation_complex_product_conic(m::JuMP.Model, a::JuMP.VariableRef, b: end - +"constraint: `c^2 + d^2 <= a*b`" function relaxation_complex_product_conic_on_off(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef, d::JuMP.VariableRef, z::JuMP.VariableRef) a_lb, a_ub = _IM.variable_domain(a) b_lb, b_ub = _IM.variable_domain(b) @@ -25,6 +25,7 @@ function relaxation_complex_product_conic_on_off(m::JuMP.Model, a::JuMP.Variable JuMP.@constraint(m, [a/sqrt(2), b_ub/sqrt(2)*z, c, d] in JuMP.RotatedSecondOrderCone()) end +"constraint: `c^2 + d^2 <= a*b*z`" function relaxation_complex_product_on_off(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef, d::JuMP.VariableRef) #to be moved to _IM a_lb, a_ub = _IM.variable_domain(a) b_lb, b_ub = _IM.variable_domain(b) @@ -47,14 +48,14 @@ function relaxation_complex_product(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP. JuMP.@constraint(m, c^2 <= a*b) end - +"constraint: `|c| <= z*b`" function relaxation_semicont_variable_on_off(m::JuMP.Model, a::JuMP.VariableRef, z::JuMP.VariableRef) a_lb, a_ub = _IM.variable_domain(a) JuMP.@constraint(m, a <= a_ub*z) JuMP.@constraint(m, a >= a_lb*z) end - +"constraint: `|c| <= z*b`" function relaxation_variable_on_off(m::JuMP.Model, x::JuMP.VariableRef, y::JuMP.VariableRef, z::JuMP.VariableRef) x_lb, x_ub = _IM.variable_domain(x) diff --git a/src/core/variableconv.jl b/src/core/variableconv.jl index 52f917cf..1254a1bc 100644 --- a/src/core/variableconv.jl +++ b/src/core/variableconv.jl @@ -980,7 +980,7 @@ function variable_converter_internal_voltage_magnitude_sqr_ne(pm::_PM.AbstractPo end report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wconv, _PM.ids(pm, nw, :convdc_ne), wcac_ne) end - +"variable: `w_du[j]` for `j` in `candidate convdc`" function variable_voltage_slack(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=false) w_du = _PM.var(pm, nw)[:w_du] = JuMP.@variable(pm.model, [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_w_du", @@ -995,7 +995,7 @@ function variable_voltage_slack(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_defau report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wdu, _PM.ids(pm, nw, :convdc_ne), w_du) end - +"variable: cos voltage fro LPAC" function variable_cos_voltage_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) #only for lpac end diff --git a/src/core/variabledcgrid.jl b/src/core/variabledcgrid.jl index 3d9e5230..6591ec4a 100644 --- a/src/core/variabledcgrid.jl +++ b/src/core/variabledcgrid.jl @@ -1,4 +1,5 @@ +"variable: DC branch currents - not used" function variable_dcbranch_current(pm::_PM.AbstractPowerModel; kwargs...) end @@ -19,7 +20,7 @@ function variable_dcgrid_voltage_magnitude(pm::_PM.AbstractPowerModel; nw::Int=_ report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :busdc, :vm, _PM.ids(pm, nw, :busdc), vdcm) end - +"variable: `vdcm[i]` for `i` in `dcbus`es" function variable_dcgrid_voltage_magnitude(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true, report::Bool=true) phivdcm = _PM.var(pm, nw)[:phi_vdcm] = JuMP.JuMP.@variable(pm.model, [i in _PM.ids(pm, nw, :busdc)], base_name="$(nw)_phi_vdcm", @@ -99,11 +100,11 @@ end ####################### TNEP variables ################################ - +"variable: DC branch currents - not used" function variable_dcbranch_current_ne(pm::_PM.AbstractPowerModel; kwargs...) end -"variable: `vdcm[i]` for `i` in `dcbus`es" +"variable: `vdcm[i]` for `i` in `ne_dcbus`es" function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) vdcm_ne = _PM.var(pm, nw)[:vdcm_ne] = JuMP.@variable(pm.model, [i in _PM.ids(pm, nw, :busdc_ne)], base_name="$(nw)_vdcm_ne", @@ -118,7 +119,7 @@ function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractPowerModel; nw::In report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :busdc_ne, :vm, _PM.ids(pm, nw, :busdc_ne), vdcm_ne) end - +"variable: `vdcm[i]` for `i` in `ne_dcbus`es" function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) phivdcm_ne = _PM.var(pm, nw)[:phi_vdcm_ne] = JuMP.@variable(pm.model, [i in _PM.ids(pm, nw, :busdc_ne)], base_name="$(nw)_phi_vdcm_ne", @@ -163,7 +164,7 @@ function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractLPACModel; nw::Int report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :branchdc_ne, :phivdcm_to, _PM.ids(pm, nw, :branchdc_ne), phivdcm_to_ne) end -"variable: `vdcm[i]` for `i` in `dcbus`es" +"variable: `wdcm[i]` for `i` in `ne_dcbus`es" function variable_dcgrid_voltage_magnitude_sqr_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) bi_bp = Dict([(i, (b["fbusdc"], b["tbusdc"])) for (i,b) in _PM.ref(pm, nw, :branchdc_ne)]) bus_vdcmax = merge(Dict([(b,bus["Vdcmax"]) for (b,bus) in _PM.ref(pm, nw, :busdc)]), @@ -192,7 +193,7 @@ function variable_dcgrid_voltage_magnitude_sqr_ne(pm::_PM.AbstractPowerModel; nw report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :busdc_ne, :wdc_ne, _PM.ids(pm, nw, :busdc_ne), wdc_ne) report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :branchdc_ne, :wdcr_ne, _PM.ids(pm, nw, :branchdc_ne), wdcr_ne) end - +"variable: `wdcm[i]` for `i` in `ne_dcbus`es" function variable_dcgrid_voltage_magnitude_sqr_du(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) # this has to to every branch, different than its counterpart(Wdc_fr) since two candidate branches can be connected to same node and two duplicate variables will be needed bi_bp = Dict([(i, (b["fbusdc"], b["tbusdc"])) for (i,b) in _PM.ref(pm, nw, :branchdc_ne)]) wdc_fr_ne = _PM.var(pm, nw)[:wdc_du_fr] = JuMP.@variable(pm.model, diff --git a/src/prob/sacdcpf.jl b/src/prob/sacdcpf.jl new file mode 100644 index 00000000..10fc738e --- /dev/null +++ b/src/prob/sacdcpf.jl @@ -0,0 +1,953 @@ +export run_sacdcpf + +""" +internal data required used solving a dc power flow + +the primary use of this data structure is to prevent re-allocation of memory +between successive power flow solves + +* `data` -- a power models data dictionary +* `busdc_gens` -- for each busdc id, a list of active generators representing dc side injections +* `amdc` -- an admittance matrix computed from the data dictionary for dc grid +* `busdc_type_idx` -- busdc types (i.e., 1, 2, 3) assigned for dc grid power flow +* `pdc_delta_base_idx` -- fixed active power delta at a busdc +* `pdc_inject_idx` -- variable active power generator injection at a busdc +* `vmdc_idx` -- variable voltage magnitude at a busdc +* `neighbors` -- neighboring buses to a given busdc +* `x0` -- 1*|N| variables, one for each busdc, varies based on busdc type # TO DO check order for dc grid +* `F0` -- 1*|N| busdc power balance evaluation values, active power only # TO DO check order for dc grid +* `J0` -- a sparse matrix holding the Jacobian of the F0 power balance evaluation function + +The postfix `_idx` indicates the admittance matrix indexing convention. +""" +struct DCPowerFlowData + data::Dict{String,<:Any} + busdc_gens::Dict{Int,Vector} + amdc::_PM.AdmittanceMatrix{Float64} + busdc_type_idx::Vector{Int} + pdc_delta_base_idx::Vector{Float64} + pdc_inject_idx::Vector{Float64} + vmdc_idx::Vector{Float64} + neighbors::Vector{Set{Int}} + x0::Vector{Float64} + F0::Vector{Float64} + J0::SparseArrays.SparseMatrixCSC{Float64,Int} +end + + +""" +This function solves sequential ac-dc power flow +""" +function run_sacdcpf(file::String; kwargs...) + data = _PM.parse_file(file) + PowerModelsACDC.process_additional_data!(data) + return run_sacdcpf(data::Dict{String,Any}, kwargs...) +end + + +function run_sacdcpf(data) + + data["load_ref"] = Dict{String, Any}() + data["load_ref"] = deepcopy(data["load"]) + data["gen_ref"] = Dict{String, Any}() + data["gen_ref"] = deepcopy(data["gen"]) + data["bus_ref"] = Dict{String, Any}() + data["bus_ref"] = deepcopy(data["bus"]) + + network = deepcopy(data) + + # STEP 1: Add converter injections as additional injections, e.g. generators (PV bus), or loads (PQ bus) + + add_converter_ac_injections!(data) + + # STEP 2: Calculate Initial AC power flow + result = _PM.compute_ac_pf(data) + if result["termination_status"] != true + Memento.warn(_LOGGER, "Initial ac powerflow in sequential acdc power flow does not converge.") + # exit() + end + + + conv_qnts = Dict{String, Any}() + pgrid_slacks = [] + result_sacdc_pf = Dict{String,Any}() + resultdc = Dict{String,Any}() + result_sacdc_pf["iterations"] = 0 + iteration = 1 + convergence = 1 + time_start_iteration = time() + while convergence > 0 + + # STEP 3: Calculate converter voltages, currents, losses, and DC side PowerModels + + conv_qnts = calc_converter_quantities(result["solution"], data) + + + vm_p = Float64[] + va_p = Float64[] + for (i,bus) in data["bus"] + push!(vm_p, result["solution"]["bus"]["$i"]["vm"]) + push!(va_p, result["solution"]["bus"]["$i"]["va"]) + end + + + # STEP 4: Add dc converter injections as additional injections in dc grid, e.g. generators (PV bus) + + add_dc_converter_injections!(data,conv_qnts) + + # STEP 5: Calculate DC grid power flows + + try + resultdc = compute_dc_pf(data, conv_qnts) + catch exception + Memento.warn(_LOGGER, "dc powerflow in sequential acdc power flow does not converge.") + break + end + + # STEP 6: Calculate AC power injections from slack converter_quantities + + pgrid_slacks = compute_slack_converter_ac_injection(resultdc["solution"], data, conv_qnts) + + # STEP 7: Update slack converter injection + + for (gen_id, pgrid) in pgrid_slacks + data["gen"]["$gen_id"]["pg"] = pgrid + end + + # STEP 8: Re-calculate AC power flows + + try + result = _PM.compute_ac_pf(data) + catch exception + Memento.warn(_LOGGER, "ac powerflow in sequential acdc power flow does not converge.") + break + end + + # STEP 9: Check "vm" and "va" for convergence + + vm_c = Float64[] + va_c = Float64[] + for (j,bus) in data["bus"] + push!(vm_c, result["solution"]["bus"]["$j"]["vm"]) + push!(va_c, result["solution"]["bus"]["$j"]["va"]) + end + + if isapprox(vm_p,vm_c; atol = 0.00001) && isapprox(abs.(va_p), abs.(va_c); atol = 0.001) + Memento.info(_LOGGER, "Sequential acdc power flow has converged.") + break + end + + iteration += 1 + + end + + + result_sacdc_pf = generate_results(result_sacdc_pf, data, result, conv_qnts, pgrid_slacks, time_start_iteration, iteration, resultdc) + + # data cleanup ! + + # Reset load in data + + data["load"] = deepcopy(data["load_ref"]) + data["bus"] = deepcopy(data["bus_ref"]) + delete!(data, "gendc") + delete!(data, "load_ref") + delete!(data, "bus_ref") + + # delete dummy gens + + gen_ref_num = maximum([gen["index"] for (g, gen) in data["gen_ref"]]) + for (i,gen) in data["gen"] + if gen["index"] > gen_ref_num + delete!(data["gen"], "$i") + delete!(result_sacdc_pf["solution"]["gen"], "$i") + end + end + + data = deepcopy(network) + + return result_sacdc_pf + +end + + + +""" +This function adds converter injections as dummy generators and loads in the ac grid +""" +function add_converter_ac_injections!(data) + load_num = maximum([load["index"] for (l, load) in data["load"]]) + gen_num = maximum([gen["index"] for (g, gen) in data["gen"]]) + bus_load_pair = Dict(load["load_bus"] => l for (l,load) in data["load_ref"]) + load_idx = 1 + gen_idx = 1 + for (c, conv) in data["convdc"] + conv_bus = conv["busac_i"] + if conv["type_dc"] == 2 + idx = gen_num + gen_idx + data["gen"]["$idx"] = Dict{String, Any}() + data["gen"]["$idx"]["pg"] = -conv["P_g"] + data["gen"]["$idx"]["qg"] = conv["Q_g"] + data["gen"]["$idx"]["model"] = 2 + data["gen"]["$idx"]["startup"] = 0.0 + data["gen"]["$idx"]["gen_bus"] = conv_bus + data["gen"]["$idx"]["vg"] = conv["Vtar"] + data["gen"]["$idx"]["mbase"] = 100 + data["gen"]["$idx"]["index"] = idx + data["gen"]["$idx"]["cost"] = [0.0, 0.0] + data["gen"]["$idx"]["qmax"] = conv["Qacmax"] + data["gen"]["$idx"]["pmax"] = conv["Pacmax"] + data["gen"]["$idx"]["qmin"] = conv["Qacmin"] + data["gen"]["$idx"]["pmin"] = conv["Pacmin"] + data["gen"]["$idx"]["ncost"] = 2 + data["gen"]["$idx"]["type"] = "dcconv" + data["gen"]["$idx"]["gen_status"] = conv["status"] + data["gen"]["$idx"]["conv_id"] = conv["index"] + data["bus"]["$conv_bus"]["bus_type"] = 2 + gen_idx += 1 + elseif conv["type_dc"] == 1 + if haskey(bus_load_pair, conv_bus) + if conv["type_ac"] == 1 + data["load"]["$(bus_load_pair[conv_bus])"]["pd"] += -conv["P_g"] + data["load"]["$(bus_load_pair[conv_bus])"]["qd"] += -conv["Q_g"] + elseif conv["type_ac"] == 2 + # adding load for active power + data["load"]["$(bus_load_pair[conv_bus])"]["pd"] += -conv["P_g"] + data["load"]["$(bus_load_pair[conv_bus])"]["qd"] += 0.0 + # adding gen for reactive power + idx = gen_num + gen_idx + data["gen"]["$idx"] = Dict{String, Any}() + data["gen"]["$idx"]["pg"] = 0.0 + data["gen"]["$idx"]["qg"] = conv["Q_g"] + data["gen"]["$idx"]["model"] = 2 + data["gen"]["$idx"]["startup"] = 0.0 + data["gen"]["$idx"]["gen_bus"] = conv_bus + data["gen"]["$idx"]["vg"] = conv["Vtar"] + data["gen"]["$idx"]["mbase"] = 100 + data["gen"]["$idx"]["index"] = idx + data["gen"]["$idx"]["cost"] = [0.0, 0.0] + data["gen"]["$idx"]["qmax"] = conv["Qacmax"] + data["gen"]["$idx"]["pmax"] = conv["Pacmax"] + data["gen"]["$idx"]["qmin"] = conv["Qacmin"] + data["gen"]["$idx"]["pmin"] = conv["Pacmin"] + data["gen"]["$idx"]["ncost"] = 2 + data["gen"]["$idx"]["type"] = "dcconv" + data["gen"]["$idx"]["gen_status"] = conv["status"] + data["gen"]["$idx"]["conv_id"] = conv["index"] + gen_idx += 1 + end + else + data["load"]["$(load_num + load_idx)"] = Dict{String, Any}() + if conv["type_ac"] == 1 + data["load"]["$(load_num + load_idx)"]["pd"] = -conv["P_g"] + data["load"]["$(load_num + load_idx)"]["qd"] = -conv["Q_g"] + elseif conv["type_ac"] == 2 + data["load"]["$(load_num + load_idx)"]["pd"] = -conv["P_g"] + data["load"]["$(load_num + load_idx)"]["qd"] = 0.0 + # adding gen for reactive power + idx = gen_num + gen_idx + data["gen"]["$idx"] = Dict{String, Any}() + data["gen"]["$idx"]["pg"] = 0.0 + data["gen"]["$idx"]["qg"] = conv["Q_g"] + data["gen"]["$idx"]["model"] = 2 + data["gen"]["$idx"]["startup"] = 0.0 + data["gen"]["$idx"]["gen_bus"] = conv_bus + data["gen"]["$idx"]["vg"] = conv["Vtar"] + data["gen"]["$idx"]["mbase"] = 100 + data["gen"]["$idx"]["index"] = idx + data["gen"]["$idx"]["cost"] = [0.0, 0.0] + data["gen"]["$idx"]["qmax"] = conv["Qacmax"] + data["gen"]["$idx"]["pmax"] = conv["Pacmax"] + data["gen"]["$idx"]["qmin"] = conv["Qacmin"] + data["gen"]["$idx"]["pmin"] = conv["Pacmin"] + data["gen"]["$idx"]["ncost"] = 2 + data["gen"]["$idx"]["type"] = "dcconv" + data["gen"]["$idx"]["gen_status"] = conv["status"] + data["gen"]["$idx"]["conv_id"] = conv["index"] + gen_idx += 1 + end + data["load"]["$(load_num + load_idx)"]["source_id"] = Any["bus", conv_bus] + data["load"]["$(load_num + load_idx)"]["load_bus"] = conv_bus + data["load"]["$(load_num + load_idx)"]["status"] = conv["status"] + data["load"]["$(load_num + load_idx)"]["index"] = load_num + load_idx + load_idx += 1 + end + end + end +end + + +""" +This function calculates converter station power flows +""" +function calc_converter_quantities(result, data) + conv_qnts = Dict{String, Any}() + for (conv_id, conv) in data["convdc"] + # Create a dict for all converters + conv_qnts["$conv_id"] = Dict{String, Any}() + # Grid voltage + conv_bus = conv["busac_i"] + conv_qnts["$conv_id"]["vm_grid"] = result["bus"]["$conv_bus"]["vm"] + conv_qnts["$conv_id"]["va_grid"] = result["bus"]["$conv_bus"]["va"] + tm = data["convdc"]["$conv_id"]["tm"] + conv_qnts["$conv_id"]["Ugrid"] = Ugrid = (result["bus"]["$conv_bus"]["vm"]*exp(-result["bus"]["$conv_bus"]["va"]im))/tm + + # Power injections + if conv["type_dc"] == 2 + for (g, gen) in data["gen"] + if haskey(gen, "type") && gen["type"] == "dcconv" && gen["conv_id"] == conv["index"] + conv_qnts["$conv_id"]["Pgrid"] = Pgrid = result["gen"][g]["pg"] + conv_qnts["$conv_id"]["Qgrid"] = Qgrid = result["gen"][g]["qg"] + end + end + elseif conv["type_dc"] == 1 + if conv["type_ac"] == 1 + bus_load_ref_pair = Dict(load["load_bus"] => l for (l,load) in data["load_ref"]) + bus_load_pair_new = Dict(load["load_bus"] => l for (l,load) in data["load"] if !haskey(data["load_ref"],l)) + if haskey(bus_load_ref_pair, conv_bus) + conv_qnts["$conv_id"]["Pgrid"] = Pgrid = data["load"]["$(bus_load_ref_pair[conv_bus])"]["pd"] - data["load_ref"]["$(bus_load_ref_pair[conv_bus])"]["pd"] + conv_qnts["$conv_id"]["Qgrid"] = Qgrid = data["load"]["$(bus_load_ref_pair[conv_bus])"]["qd"] - data["load_ref"]["$(bus_load_ref_pair[conv_bus])"]["qd"] + elseif haskey(bus_load_pair_new, conv_bus) + conv_qnts["$conv_id"]["Pgrid"] = Pgrid = data["load"]["$(bus_load_pair_new[conv_bus])"]["pd"] + conv_qnts["$conv_id"]["Qgrid"] = Qgrid = data["load"]["$(bus_load_pair_new[conv_bus])"]["qd"] + end + elseif conv["type_ac"] == 2 + bus_load_ref_pair = Dict(load["load_bus"] => l for (l,load) in data["load_ref"]) + bus_load_pair_new = Dict(load["load_bus"] => l for (l,load) in data["load"] if !haskey(data["load_ref"],l)) + if haskey(bus_load_ref_pair, conv_bus) + conv_qnts["$conv_id"]["Pgrid"] = Pgrid = data["load"]["$(bus_load_ref_pair[conv_bus])"]["pd"] - data["load_ref"]["$(bus_load_ref_pair[conv_bus])"]["pd"] + elseif haskey(bus_load_pair_new, conv_bus) + conv_qnts["$conv_id"]["Pgrid"] = Pgrid = data["load"]["$(bus_load_pair_new[conv_bus])"]["pd"] + end + for (g, gen) in data["gen"] + if haskey(gen, "type") && gen["type"] == "dcconv" && gen["conv_id"] == conv["index"] + conv_qnts["$conv_id"]["Qgrid"] = Qgrid = result["gen"][g]["qg"] + end + end + end + + end + conv_qnts["$conv_id"]["Sgrid"] = Sgrid = Pgrid + Qgrid*im + # Transformer current calculation: Itf = Sgrid / Ugrid + Itf = conj(Sgrid / Ugrid) + # Filter current If = -Ufilter * Bf = -(Ugrid - Itf * (Rtf + j Ztf)) + Ztf = (data["convdc"]["$conv_id"]["rtf"] + data["convdc"]["$conv_id"]["xtf"]im ) * data["convdc"]["$conv_id"]["transformer"] + conv_qnts["$conv_id"]["Ztf"] = Ztf + conv_qnts["$conv_id"]["Bf"] = Bf = data["convdc"]["$conv_id"]["bf"] * data["convdc"]["$conv_id"]["filter"] + # If = -(Ugrid - Itf*(Ztf)) * Bf + conv_qnts["$conv_id"]["Uf"] = Uf = (Ugrid - Itf*(Ztf)) + conv_qnts["$conv_id"]["If"] = If = (Ugrid - Itf*(Ztf)) * (-Bf*im) + # Reactor current Ipr = Itf - If + conv_qnts["$conv_id"]["Ipr"] = Ipr = Itf - If + # Converter voltage Uc = Uf - Ic * Zpr & Ic = Ipr + Zpr = (data["convdc"]["$conv_id"]["rc"] + data["convdc"]["$conv_id"]["xc"]im ) * data["convdc"]["$conv_id"]["reactor"] + conv_qnts["$conv_id"]["Zpr"] = Zpr + Uc = (Ugrid - Itf*(Ztf)) - Ipr * Zpr + # Converter power Sconv = Uc * Ic' + conv_qnts["$conv_id"]["Sconv"] = Sconv = Uc * conj(Ipr) + Pconv = real(Sconv) + Qconv = imag(Sconv) + # Converter losses Ploss = a + b * |Ic| + c * Ic^2 + Ploss = data["convdc"]["$conv_id"]["LossA"] + data["convdc"]["$conv_id"]["LossB"] * abs(Ipr) + data["convdc"]["$conv_id"]["LossCrec"] * abs(Ipr)^2 + # Pdc = Pconv - Ploss + Pdc = -Pconv + Ploss + conv_qnts["$conv_id"]["Pconv"] = Pconv + conv_qnts["$conv_id"]["Qconv"] = Qconv + conv_qnts["$conv_id"]["Pdc"] = Pdc + conv_qnts["$conv_id"]["Ploss"] = Ploss + conv_qnts["$conv_id"]["Uc"] = Uc + conv_qnts["$conv_id"]["Ptf_to"] = real(Uf * -conj(Itf)) + conv_qnts["$conv_id"]["Qtf_to"] = imag(Uf * -conj(Itf)) + conv_qnts["$conv_id"]["Ppr_fr"] = real(Uf * conj(Ipr)) + conv_qnts["$conv_id"]["Qpr_fr"] = imag(Uf * conj(Ipr)) + + end + + return conv_qnts +end + + +""" +This function adds converter injections as dummy generators in the dc grid +""" +function add_dc_converter_injections!(data,conv_qnts) + # add dummy generators on busdc + if haskey(data, "gendc") + for (idx, conv) in data["convdc"] + if conv["type_dc"] == 2 + data["gendc"]["$idx"]["pg"] = conv_qnts["$idx"]["Pdc"] + end + end + else + data["gendc"] = Dict{String, Any}() + for (idx, conv) in data["convdc"] + data["gendc"]["$idx"] = Dict{String, Any}() + data["gendc"]["$idx"]["pg"] = -conv_qnts["$idx"]["Pdc"] + data["gendc"]["$idx"]["qg"] = 0.0 + data["gendc"]["$idx"]["model"] = 2 + data["gendc"]["$idx"]["startup"] = 0.0 + data["gendc"]["$idx"]["gen_bus"] = conv["busdc_i"] + data["gendc"]["$idx"]["vg"] = conv["Vtar"] + data["gendc"]["$idx"]["mbase"] = 100 + data["gendc"]["$idx"]["index"] = idx + data["gendc"]["$idx"]["cost"] = [0.0, 0.0] + data["gendc"]["$idx"]["qmax"] = 0.0 + data["gendc"]["$idx"]["pmax"] = 1.2*conv["Pacmax"] + data["gendc"]["$idx"]["qmin"] = 0.0 + data["gendc"]["$idx"]["pmin"] = conv["Pacmin"] + data["gendc"]["$idx"]["ncost"] = 2 + data["gendc"]["$idx"]["type"] = "dcconv_busdc" + data["gendc"]["$idx"]["gen_status"] = conv["status"] + data["gendc"]["$idx"]["conv_id"] = conv["index"] + bus_idx = conv["busdc_i"] + if conv["type_dc"] == 1 || conv["type_dc"] == 3 + data["busdc"]["$bus_idx"]["bus_type"] = 1 + elseif conv["type_dc"] == 2 + data["busdc"]["$bus_idx"]["bus_type"] = 3 + end + end + end +end + + +""" +This function solves dc grid power flow +""" +function compute_dc_pf(data, conv_qnts) + + dcpf_data = instantiate_dcpf_data(data, conv_qnts) + dc_pf_results = _compute_dc_pf(dcpf_data) + + return dc_pf_results + +end + + +function instantiate_dcpf_data(data::Dict{String,<:Any}, conv_qnts::Dict{String,<:Any}) + + pdc_delta = calc_busdc_injection(data) + + # remove gendc injections from slack + for (i,gendc) in data["gendc"] + gendc_bus = data["busdc"]["$(gendc["gen_bus"])"] + if gendc["gen_status"] != 0 + if gendc_bus["bus_type"] == 3 + pdc_delta[gendc_bus["index"]] -= gendc["pg"] + end + end + end + + + busdc_gens = Dict{Int,Array{Any}}() + for (i,gendc) in data["gendc"] + # skip inactive generators + if gendc["gen_status"] == 0 + continue + end + + gendc_bus_id = gendc["gen_bus"] + if !haskey(busdc_gens, gendc_bus_id) + busdc_gens[gendc_bus_id] = [] + end + push!(busdc_gens[gendc_bus_id], gendc) + end + + for (busdc_id, gensdc) in busdc_gens + sort!(gensdc, by=x -> (x["qmax"] - x["qmin"], x["index"])) + end + + amdc = calc_admittance_matrix(data) + + busdc_type_idx = Int[data["busdc"]["$(bus_id)"]["bus_type"] for bus_id in amdc.idx_to_bus] + + pdc_delta_base_idx = Float64[-pdc_delta[bus_id] for bus_id in amdc.idx_to_bus] + + pdc_inject_idx = [0.0 for bus_id in amdc.idx_to_bus] + + vmdc_idx = [1.0 for bus_id in amdc.idx_to_bus] + + + # for buses with non-1.0 bus voltages + for (i,busdc) in data["busdc"] + if busdc["bus_type"] == 2 || busdc["bus_type"] == 3 + vmdc_idx[amdc.bus_to_idx[busdc["index"]]] = busdc["Vdc"] + end + end + + + neighbors = [Set{Int}([i]) for i in eachindex(amdc.idx_to_bus)] + I, J, V = _PM.findnz(amdc.matrix) + for nz in eachindex(V) + push!(neighbors[I[nz]], J[nz]) + push!(neighbors[J[nz]], I[nz]) + end + + x0 = [0.0 for i in 1:length(amdc.idx_to_bus)] + F0 = similar(x0) + + J0_I = Int[] + J0_J = Int[] + J0_V = Float64[] + + for i in eachindex(amdc.idx_to_bus) + + for j in neighbors[i] + push!(J0_I, i); push!(J0_J, j); push!(J0_V, 0.0) + end + end + J0 = _PM.sparse(J0_I, J0_J, J0_V) + + return DCPowerFlowData(data, busdc_gens, amdc, busdc_type_idx, pdc_delta_base_idx, pdc_inject_idx, vmdc_idx, neighbors, x0, F0, J0) +end + + +function _compute_dc_pf(dcpf_data::DCPowerFlowData; finite_differencing=false, flat_start=false, kwargs...) + time_start = time() + data = dcpf_data.data + amdc = dcpf_data.amdc + busdc_type_idx = dcpf_data.busdc_type_idx + pdc_delta_base_idx = dcpf_data.pdc_delta_base_idx + pdc_inject_idx = dcpf_data.pdc_inject_idx + vmdc_idx = dcpf_data.vmdc_idx + neighbors = dcpf_data.neighbors + x0 = dcpf_data.x0 + F0 = dcpf_data.F0 + J0 = dcpf_data.J0 + + # dc power flow, nodal power balance function eval + function f!(F::Vector{Float64}, x::Vector{Float64}) + for i in eachindex(amdc.idx_to_bus) + if busdc_type_idx[i] == 1 + vmdc_idx[i] = x[i] + elseif busdc_type_idx[i] == 2 + elseif busdc_type_idx[i] == 3 + pdc_inject_idx[i] = x[i] + else + @assert false + end + end + + for i in eachindex(amdc.idx_to_bus) + balance_real = pdc_delta_base_idx[i] + pdc_inject_idx[i] + for j in neighbors[i] + if i == j + balance_real += vmdc_idx[i] * vmdc_idx[i] * amdc.matrix[i,i] + else + balance_real += vmdc_idx[i] * vmdc_idx[j] * amdc.matrix[i,j] + end + end + F[i] = balance_real + end + + + end + + + # dc power flow, sparse jacobian computation + function jsp!(J::_PM.SparseMatrixCSC{Float64,Int}, x::Vector{Float64}) + for i in eachindex(amdc.idx_to_bus) + + for j in neighbors[i] + + bus_type = busdc_type_idx[j] + if bus_type == 1 + if i == j + y_ii = amdc.matrix[i,i] + J[i, j] = 2 * y_ii * vmdc_idx[i] + sum(amdc.matrix[i,k] * vmdc_idx[k] for k in neighbors[i] if k != i) + else + y_ij = amdc.matrix[i,j] + J[i, j] = vmdc_idx[i] * y_ij + end + elseif bus_type == 2 + if i == j + J[i, j] = 1.0 + else + J[i, j] = 0.0 + end + elseif bus_type == 3 + if i == j + J[i, j] = 1.0 + else + J[i, j] = 0.0 + end + else + @assert false + end + end + end + end + + + # basic init point + for i in eachindex(amdc.idx_to_bus) + if busdc_type_idx[i] == 1 + x0[i] = 1.0 + elseif busdc_type_idx[i] == 2 + elseif busdc_type_idx[i] == 3 + else + @assert false + end + end + + # warm-start point + if !flat_start + pdc_inject = Dict{Int,Float64}(busdc["index"] => 0.0 for (i,busdc) in data["busdc"]) + for (i,gendc) in data["gendc"] + if gendc["gen_status"] != 0 + if haskey(gendc, "pg_start") + pdc_inject[gendc["gen_bus"]] += gendc["pg_start"] + end + end + end + + + for (i,bid) in enumerate(amdc.idx_to_bus) + busdc = data["busdc"]["$(bid)"] + if busdc_type_idx[i] == 1 + if haskey(busdc, "vm_start") + x0[i] = busdc["vm_start"] + end + elseif busdc_type_idx[i] == 2 + elseif busdc_type_idx[i] == 3 + x0[i] = -pdc_inject[bid] + else + @assert false + end + end + end + + + + if finite_differencing + result = NLsolve.nlsolve(f!, x0; kwargs...) + else + df = NLsolve.OnceDifferentiable(f!, jsp!, x0, F0, J0) + result = NLsolve.nlsolve(df, x0, xtol = 0.0, ftol = 1e-8, iterations = 1000, show_trace = false) + end + + solution = Dict("per_unit" => dcpf_data.data["per_unit"]) + + converged = result.x_converged || result.f_converged + + if !converged + Memento.warn(_LOGGER, "dc power flow solver convergence failed! use `show_trace = true` for more details") + else + data = dcpf_data.data + busdc_gens = dcpf_data.busdc_gens + amdc = dcpf_data.amdc + busdc_type_idx = dcpf_data.busdc_type_idx + + busdc_assignment= Dict{String,Any}() + for (i,busdc) in data["busdc"] + if busdc["bus_type"] != 4 + busdc_idx = amdc.bus_to_idx[busdc["index"]] + + busdc_assignment[i] = Dict( + "Vdc" => dcpf_data.vmdc_idx[busdc_idx] + ) + end + end + + gendc_assignment= Dict{String,Any}() + for (i,gendc) in data["gendc"] + if gendc["gen_status"] != 0 + gendc_assignment[i] = Dict( + "pg" => gendc["pg"] + ) + end + end + + + for (i,bid) in enumerate(amdc.idx_to_bus) + busdc = busdc_assignment["$(bid)"] + + if busdc_type_idx[i] == 1 + busdc["Vdc"] = result.zero[i] + elseif busdc_type_idx[i] == 2 + for gendc in busdc_gens[bid] + sol_gen = gendc_assignment["$(gendc["index"])"] + end + elseif busdc_type_idx[i] == 3 + for gendc in busdc_gens[bid] + sol_gen = gendc_assignment["$(gendc["index"])"] + sol_gen["pg"] = 0.0 + end + pg_remaining = result.zero[i] + _assign_pg!(gendc_assignment, busdc_gens[bid], pg_remaining) + else + @assert false + end + end + + solution = Dict( + "per_unit" => data["per_unit"], + "busdc" => busdc_assignment, + "gendc" => gendc_assignment, + ) + end + + results = Dict( + "optimizer" => "NLsolve", + "termination_status" => converged, + "objective" => 0.0, + "solution" => solution, + "solve_time" => time() - time_start + ) + + return results +end + + +""" +This function calculates dc bus injections +""" +function calc_busdc_injection(data::Dict{String,<:Any}) + busdc_values = Dict(busdc["index"] => Dict{String,Float64}() for (i,busdc) in data["busdc"]) + for (i,busdc) in data["busdc"] + bvals = busdc_values[busdc["index"]] + bvals["pg"] = 0.0 + end + for (i,gendc) in data["gendc"] + if gendc["gen_status"] != 0 + bvals = busdc_values[gendc["gen_bus"]] + bvals["pg"] += gendc["pg"] + end + end + pdc_delta = Dict{Int,Float64}() + for (i,busdc) in data["busdc"] + if busdc["bus_type"] != 4 + bvals = busdc_values[busdc["index"]] + p_delta = bvals["pg"] + else + p_delta = NaN + end + pdc_delta[busdc["index"]] = p_delta + end + return pdc_delta +end + + +""" +This function calculates dc admittance matrix +""" +function calc_admittance_matrix(data::Dict{String,<:Any}) + + dc_buses = [x.second for x in data["busdc"]] + sort!(dc_buses, by=x->x["index"]) + + idx_to_busdc = [x["index"] for x in dc_buses] + busdc_to_idx = Dict(x["index"] => i for (i,x) in enumerate(dc_buses)) + + I = Int[] + J = Int[] + V = Float64[] + + for (i,branchdc) in data["branchdc"] + f_bus = branchdc["fbusdc"] + t_bus = branchdc["tbusdc"] + if branchdc["status"] != 0 && haskey(busdc_to_idx, f_bus) && haskey(busdc_to_idx, t_bus) + f_bus = busdc_to_idx[f_bus] + t_bus = busdc_to_idx[t_bus] + g = inv(branchdc["r"]) + p = data["dcpol"] + push!(I, f_bus); push!(J, t_bus); push!(V, -p*g) + push!(I, t_bus); push!(J, f_bus); push!(V, -p*g) + push!(I, f_bus); push!(J, f_bus); push!(V, p*g) + push!(I, t_bus); push!(J, t_bus); push!(V, p*g) + end + end + + m = _PM.sparse(I,J,V) + + amdc = _PM.AdmittanceMatrix(idx_to_busdc, busdc_to_idx, m) + return amdc +end + + +""" +This function performs internal iteration to calculate slack converter ac grid active injection +""" +function compute_slack_converter_ac_injection(resultdc, data, conv_qnts) + Pdc1 = 0.0 + pgrid_slacks_new = [] + slack_conv_busac_i = 0 + for (conv_id, conv) in data["convdc"] + if conv["type_dc"] == 2 + for (gen_id, gendc) in data["gendc"] + if gendc["gen_bus"] == conv["busdc_i"] + Pdc1 = -resultdc["gendc"]["$gen_id"]["pg"] + end + end + + Ploss0 = conv_qnts["$conv_id"]["Ploss"] + Qconv0 = conv_qnts["$conv_id"]["Qconv"] + Uc0 = conv_qnts["$conv_id"]["Uc"] + Ugrid0 = conv_qnts["$conv_id"]["Ugrid"] + Zpr = conv_qnts["$conv_id"]["Zpr"] + Ztf = conv_qnts["$conv_id"]["Ztf"] + Bf = conv_qnts["$conv_id"]["Bf"] + + # iteration 1-Fwd + Pconv1 = -Pdc1 + Ploss0 + Sconv1 = Pconv1 + Qconv0*im + Ipr1 = conj(Sconv1/Uc0) + If1 = (Uc0 + Ipr1 * Zpr) * (-Bf*im) + Itf1 = If1 + Ipr1 + Sgrid1 = Ugrid0 * conj(Itf1) + + # iteration 1-Rev + Itf2 = conj(Sgrid1/Ugrid0) + # If1 = -(Ugrid0 - Itf2*(Ztf)) * Bf + If2 = (Ugrid0 - Itf2*(Ztf)) * (-Bf*im) + Ipr2 = Itf2 - If2 + # Uc1 = (Ugrid0 - Itf2*(Ztf)) - Ipr2 * Zpr + Sconv2 = Uc0 * conj(Ipr2) + Pconv2 = real(Sconv2) + Ploss1 = Pconv2 + Pdc1 + + # Extract new values + Pconv_new = -Pdc1 + Ploss1 + Sconv_new = Pconv_new + Qconv0*im + Ipr_new = conj(Sconv_new/Uc0) + If_new = (Uc0 + Ipr2 * Zpr) * (-Bf*im) + Itf_new = If_new + Ipr_new + Sgrid_new = Ugrid0 * conj(Itf_new) + slack_conv_busac_i = conv["busac_i"] + + # + conv_qnts["$conv_id"]["Pgrid"] = real(Sgrid_new) + conv_qnts["$conv_id"]["Qgrid"] = imag(Sgrid_new) + + for (g, gen) in data["gen"] + if haskey(gen, "type") && gen["type"] == "dcconv" && gen["gen_bus"] == slack_conv_busac_i && gen["conv_id"] == conv["index"] + push!( pgrid_slacks_new, (gen["index"], real(Sgrid_new)) ) + end + end + + end + end + + return pgrid_slacks_new +end + + +""" +This function generates results similar to acdcpf +""" +function generate_results(result_sacdc_pf, data, result, conv_qnts, pgrid_slacks, time_start_iteration, iteration, resultdc) + result_sacdc_pf["time_iteration"] = time() - time_start_iteration + result_sacdc_pf["iterations"] = iteration + result_sacdc_pf["solution"] = result["solution"] + result_sacdc_pf["optimizer"] = NLsolve + result_sacdc_pf["per_unit"] = data["per_unit"] + result_sacdc_pf["termination_status"] = "Converged" + result_sacdc_pf["objective"] = 0.0 + result_sacdc_pf["solution"]["convdc"] = Dict{String,Any}() + + for (i,conv) in data["convdc"] + if conv["status"] != 0 + if conv["type_dc"] != 2 + result_sacdc_pf["solution"]["convdc"][i] = Dict( + "vmfilt" => abs(conv_qnts[i]["Uf"]), + "qpr_fr" => conv_qnts[i]["Qpr_fr"], + "ppr_fr" => conv_qnts[i]["Ppr_fr"], + "qconv" => imag(conv_qnts[i]["Sconv"]), + "iconv" => abs(conv_qnts[i]["Ipr"]), + "pgrid" => -data["convdc"][i]["P_g"], + "qtf_to" => conv_qnts[i]["Qtf_to"], + "phi" => angle(conv_qnts[i]["Sconv"]), + "vaconv" => angle(conv_qnts[i]["Uc"]), + "pconv" => conv_qnts[i]["Pconv"], + "ptf_to" => conv_qnts[i]["Ptf_to"], + "vmconv" => abs(conv_qnts[i]["Uc"]), + "vafilt" => angle(conv_qnts[i]["Uf"]), + "pdc" => conv_qnts[i]["Pdc"], + "qgrid" => -data["convdc"][i]["Q_g"] + ) + bus_load_ref_pair = Dict(load["load_bus"] => l for (l,load) in data["load_ref"]) + bus_load_pair_new = Dict(load["load_bus"] => l for (l,load) in data["load"] if !haskey(data["load_ref"],l)) + conv_bus = conv["busac_i"] + if haskey(bus_load_ref_pair, conv_bus) + result_sacdc_pf["solution"]["convdc"][i]["pgrid"] = - (data["load"]["$(bus_load_ref_pair[conv_bus])"]["pd"] - data["load_ref"]["$(bus_load_ref_pair[conv_bus])"]["pd"]) + result_sacdc_pf["solution"]["convdc"][i]["qgrid"] = - (data["load"]["$(bus_load_ref_pair[conv_bus])"]["qd"] - data["load_ref"]["$(bus_load_ref_pair[conv_bus])"]["qd"]) + elseif haskey(bus_load_pair_new, conv_bus) + result_sacdc_pf["solution"]["convdc"][i]["pgrid"] = - (data["load"]["$(bus_load_pair_new[conv_bus])"]["pd"]) + result_sacdc_pf["solution"]["convdc"][i]["qgrid"] = - (data["load"]["$(bus_load_pair_new[conv_bus])"]["qd"]) + end + # + if conv["type_ac"] == 2 + result_sacdc_pf["solution"]["convdc"][i]["qgrid"] = [result["solution"]["gen"][g]["qg"] for (g, gen) in data["gen"] if haskey(gen, "type") && gen["type"] == "dcconv" && gen["conv_id"] == conv["index"]][1] + end + else + result_sacdc_pf["solution"]["convdc"][i] = Dict( + "vmfilt" => abs(conv_qnts[i]["Uf"]), + "qpr_fr" => conv_qnts[i]["Qpr_fr"], + "ppr_fr" => conv_qnts[i]["Ppr_fr"], + "qconv" => imag(conv_qnts[i]["Sconv"]), + "iconv" => abs(conv_qnts[i]["Ipr"]), + "pgrid" => [result["solution"]["gen"][g]["pg"] for (g, gen) in data["gen"] if haskey(gen, "type") && gen["type"] == "dcconv" && gen["conv_id"] == conv["index"]][1], + "qtf_to" => conv_qnts[i]["Qtf_to"], + "phi" => angle(conv_qnts[i]["Sconv"]), + "vaconv" => angle(conv_qnts[i]["Uc"]), + "pconv" => real(conv_qnts[i]["Sconv"]), + "ptf_to" => conv_qnts[i]["Ptf_to"], + "vmconv" => abs(conv_qnts[i]["Uc"]), + "vafilt" => angle(conv_qnts[i]["Uf"]), + "pdc" => -conv_qnts[i]["Pdc"], + "qgrid" => [result["solution"]["gen"][g]["qg"] for (g, gen) in data["gen"] if haskey(gen, "type") && gen["type"] == "dcconv" && gen["conv_id"] == conv["index"]][1] + ) + end + end + end + + result_sacdc_pf["solution"]["busdc"] = Dict{String,Any}() + for (i,bsdc) in data["busdc"] + result_sacdc_pf["solution"]["busdc"][i] = Dict( + "vm" => resultdc["solution"]["busdc"][i]["Vdc"] + ) + end + + result_sacdc_pf["solution"]["branchdc"] = Dict{String,Any}() + for (i,brdc) in data["branchdc"] + if brdc["status"] != 0 + result_sacdc_pf["solution"]["branchdc"][i] = Dict( + "pt" => data["dcpol"] * (1/brdc["r"]) * resultdc["solution"]["busdc"]["$(brdc["tbusdc"])"]["Vdc"] * (resultdc["solution"]["busdc"]["$(brdc["tbusdc"])"]["Vdc"] - resultdc["solution"]["busdc"]["$(brdc["fbusdc"])"]["Vdc"]), + "pf" => data["dcpol"] * (1/brdc["r"]) * resultdc["solution"]["busdc"]["$(brdc["fbusdc"])"]["Vdc"] * (resultdc["solution"]["busdc"]["$(brdc["fbusdc"])"]["Vdc"] - resultdc["solution"]["busdc"]["$(brdc["tbusdc"])"]["Vdc"]) + ) + end + end + return result_sacdc_pf +end + + +""" +This function assigns active power to active generators representing dc side injections within limits +""" +function _assign_pg!(sol_gens::Dict{String,<:Any}, busdc_gens::Vector, pg_remaining::Float64) + for gendc in busdc_gens[1:end-1] + pmin = gendc["pmin"] + pmax = gendc["pmax"] + + if (pg_remaining <= 0.0 && pmin >= 0.0) || (pg_remaining >= 0.0 && pmax <= 0.0) + # keep pg assignment as zero + continue + end + + sol_gen = sol_gens["$(gendc["index"])"] + if pg_remaining < pmin + sol_gen["pg"] = pmin + elseif pg_remaining > pmax + sol_gen["pg"] = pmax + else + sol_gen["pg"] = pg_remaining + pg_remaining = 0.0 + break + end + pg_remaining -= sol_gen["pg"] + end + if !isapprox(pg_remaining, 0.0) + gendc = busdc_gens[end] + sol_gen = sol_gens["$(gendc["index"])"] + sol_gen["pg"] = pg_remaining + end +end diff --git a/test/runtests.jl b/test/runtests.jl index c25783fc..dacfb97c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,4 +44,6 @@ include("opf.jl") include("tnep.jl") +include("spf.jl") + end diff --git a/test/scripts/ACDCtest.jl b/test/scripts/ACDCtest.jl index 69837e5e..d9be630a 100644 --- a/test/scripts/ACDCtest.jl +++ b/test/scripts/ACDCtest.jl @@ -31,7 +31,7 @@ s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true) result = _PMACDC.run_acdcpf(file, _PM.ACRPowerModel, ipopt; setting = s) result_droop = _PMACDC.run_acdcpf(file_pf_droop, _PM.ACRPowerModel, ipopt; setting = s) - +result = _PMACDC.run_acdcpf(file_case3120, _PM.ACRPowerModel, ipopt; setting = s) # resultAC = _PMACDC.run_acdcopf(file, _PM.ACPPowerModel, ipopt; setting = s) # resultLPAC = _PMACDC.run_acdcopf(file, _PM.LPACCPowerModel, ipopt; setting = s) diff --git a/test/spf.jl b/test/spf.jl new file mode 100644 index 00000000..eb03a0ec --- /dev/null +++ b/test/spf.jl @@ -0,0 +1,26 @@ +@testset "test sacdc pf" begin + @testset "5-bus ac dc case" begin + result = run_sacdcpf("../test/data/case5_acdc.m") + + @test result["termination_status"] == "Converged" + @test isapprox(result["objective"], 0.0; atol = 1e-2) + + @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.32323; atol = 2e-3) + @test isapprox(result["solution"]["bus"]["5"]["va"], -0.07172; atol = 2e-3) + @test isapprox(result["solution"]["bus"]["3"]["va"], -0.06676; atol = 2e-3) + @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], 0.22012; atol = 2e-3) + @test isapprox(result["solution"]["convdc"]["3"]["pdc"], 0.36422; atol = 2e-3) + end + @testset "39-bus ac dc case" begin + result = run_sacdcpf("../test/data/case39_acdc.m") + + @test result["termination_status"] == "Converged" + @test isapprox(result["objective"], 0.0; atol = 1e-2) + + @test isapprox(result["solution"]["gen"]["7"]["pg"], 5.60000; atol = 2e-3) + @test isapprox(result["solution"]["bus"]["20"]["va"], -0.39552; atol = 2e-3) + @test isapprox(result["solution"]["bus"]["27"]["va"], -0.48778; atol = 2e-3) + @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.60000; atol = 2e-3) + @test isapprox(result["solution"]["convdc"]["3"]["pdc"], -0.58734; atol = 2e-3) + end +end \ No newline at end of file