added pflex

Project.toml

@@ -1,7 +1,7 @@
 name = "PandaModels"
 uuid = "2dbab86a-7cbf-476f-9afe-75ffd3079e7c"
 authors = ["e2nIEE"]
-version = "0.7.2"
+version = "0.7.3"
 
 [deps]
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"

src/PandaModels.jl

@@ -28,6 +28,7 @@ export run_powermodels_pf,
     run_powermodels_multi_storage,
     run_pandamodels_multi_vstab,
     run_pandamodels_qflex,
+    run_pandamodels_pflex,
     run_pandamodels_multi_qflex,
     run_pandamodels_ploss,
     run_pandamodels_vstab_test,
@@ -39,6 +40,7 @@ include("input/pp_to_pm.jl")
 include("input/tools.jl")
 include("models/vstab.jl")
 include("models/qflex.jl")
+include("models/pflex.jl")
 include("models/ploss.jl")
 include("models/call_pandamodels.jl")
 include("models/call_powermodels.jl")

src/models/call_pandamodels.jl

@@ -50,6 +50,24 @@ function run_pandamodels_qflex(json_path)
     return result
 end
 
+function run_pandamodels_pflex(json_path)
+    pm = load_pm_from_json(json_path)
+    active_powermodels_silence!(pm)
+    pm = check_powermodels_data!(pm)
+    model = get_model(pm["pm_model"])
+    solver = get_solver(pm)
+
+    result = _run_pflex(
+        pm,
+        model,
+        solver,
+        setting = Dict("output" => Dict("branch_flows" => true)),
+        ext = extract_params!(pm),
+    )
+    return result
+end
+
+
 function run_pandamodels_multi_qflex(json_path)
     pm = load_pm_from_json(json_path)
     active_powermodels_silence!(pm)

src/models/pflex.jl

@@ -0,0 +1,118 @@
+export _run_pflex, _run_multi_pflex
+
+"""
+run optimization for maintaining power setpoints
+"""
+
+function _run_pflex(file, model_type::_PM.Type, optimizer; kwargs...)
+    return _PM.solve_model(file, model_type, optimizer, _build_pflex; kwargs...)
+end
+
+"""
+give a JuMP model with PowerModels network data structur and build opitmization model
+"""
+
+function _build_pflex(pm::_PM.AbstractPowerModel)
+
+    _PM.variable_bus_voltage(pm)
+    _PM.variable_gen_power(pm)
+    _PM.variable_branch_power(pm)
+    _PM.variable_dcline_power(pm, bounded = false) # TODO: why false?
+
+    objective_pflex(pm)
+    # println("pflex objective function:", JuMP.objective_function(pm.model))
+    _PM.constraint_model_voltage(pm)
+
+    for i in _PM.ids(pm, :ref_buses)
+        _PM.constraint_theta_ref(pm, i)
+    end
+
+    for i in _PM.ids(pm, :bus)
+        _PM.constraint_power_balance(pm, i)
+    end
+
+    for (i, branch) in _PM.ref(pm, :branch)
+        _PM.constraint_ohms_yt_from(pm, i)
+        _PM.constraint_ohms_yt_to(pm, i)
+        _PM.constraint_thermal_limit_from(pm, i)
+        _PM.constraint_thermal_limit_to(pm, i)
+    end
+
+    for i in _PM.ids(pm, :dcline)
+        _PM.constraint_dcline_power_losses(pm, i)
+    end
+end
+
+function objective_pflex(pm::_PM.AbstractPowerModel)
+
+    if haskey(pm.ext, :obj_factors)
+        if length(pm.ext[:obj_factors]) == 2
+            fac1 = pm.ext[:obj_factors]["fac_1"]
+            fac2 = pm.ext[:obj_factors]["fac_2"]
+        end
+    else
+        fac1 = 1.0
+        fac2 = 0
+    end
+
+    return JuMP.@objective(pm.model, Min,
+    fac1 * sum((var(pm, :p, (content["element_index"], content["f_bus"], content["t_bus"])) - content["value"])^2
+        for (i, content) in pm.ext[:setpoint_p])
+    +
+    fac2 * sum((var(pm, :pg, content)-0)^2 for (i, content) in pm.ext[:gen_and_controllable_sgen]))
+end
+
+
+function _run_multi_pflex(file, model_type::_PM.Type, optimizer; kwargs...)
+    return _PM.solve_model(file, model_type, optimizer, _build_multi_pflex; multinetwork=true, kwargs...)
+end
+
+"""
+run multi-timestep optimization for maintaining power setpoints
+"""
+
+function _build_multi_pflex(pm::_PM.AbstractPowerModel)
+    for (n, network) in _PM.nws(pm)
+            _PM.variable_bus_voltage(pm, nw=n)
+            _PM.variable_gen_power(pm, nw=n)
+            _PM.variable_branch_power(pm, nw=n)
+            _PM.variable_dcline_power(pm, nw=n)
+
+            _PM.constraint_model_voltage(pm, nw=n)
+
+            for i in ids(pm, :ref_buses, nw=n)
+                _PM.constraint_theta_ref(pm, i, nw=n)
+            end
+
+            for i in ids(pm, :bus, nw=n)
+                _PM.constraint_power_balance(pm, i, nw=n)
+            end
+
+            for i in ids(pm, :branch, nw=n)
+                _PM.constraint_ohms_yt_from(pm, i, nw=n)
+                _PM.constraint_ohms_yt_to(pm, i, nw=n)
+
+                _PM.constraint_voltage_angle_difference(pm, i, nw=n)
+
+                _PM.constraint_thermal_limit_from(pm, i, nw=n)
+                _PM.constraint_thermal_limit_to(pm, i, nw=n)
+            end
+
+            for i in ids(pm, :dcline, nw=n)
+                _PM.constraint_dcline_power_losses(pm, i, nw=n)
+            end
+    end
+
+    objective_multi_pflex(pm)
+
+end
+
+
+function objective_multi_pflex(pm::_PM.AbstractPowerModel)
+    timestep_ids = [id for id in _PM.nw_ids(pm) if id != 0]
+    return JuMP.@objective(pm.model, Min,
+        sum(
+        sum((var(pm, nw, :p, (content["element_index"], content["f_bus"], content["t_bus"])) - content["value"])^2 for (i, content) in pm.ext[:setpoint_p])
+        for nw in timestep_ids)
+            )
+end

test/call_pandamodels.jl

@@ -39,6 +39,24 @@
         @test result["solve_time"] > 0.0
     end
 
+    @testset "case_pflex: cigre mv" begin
+        result = run_pandamodels_pflex(case_pflex)
+        pm = _PdM.load_pm_from_json(case_pflex)
+        params = _PdM.extract_params!(pm)
+        print(params[:setpoint_p])
+        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
+        @test string(result["dual_status"]) == "FEASIBLE_POINT"
+        @test string(result["primal_status"]) == "FEASIBLE_POINT"
+
+        for (idx, br) in result["solution"]["branch"]
+            if idx in keys(params[:setpoint_p])
+                @test isapprox(br["pf"], params[:setpoint_p][idx]["value"], atol=1e-1)
+            end
+        end
+        @test isapprox(result["objective_lb"], -Inf)
+        @test result["solve_time"] > 0.0
+    end
+
     @testset "case_multi_qflex: cigre mv" begin
         result = run_pandamodels_multi_qflex(case_multi_qflex)
         pm = _PdM.load_pm_from_json(case_multi_qflex)