|
| 1 | +function run_pandamodels_qflex_test(json_path) # before run_poweramodels |
| 2 | + time_start = time() |
| 3 | + ############################################################################### |
| 4 | + # 0. Initialization |
| 5 | + ############################################################################### |
| 6 | + pm = load_pm_from_json(json_path) |
| 7 | + solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"], |
| 8 | + pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"]) |
| 9 | + if haskey(pm, "user_defined_params") |
| 10 | + user_defined_params = pm["user_defined_params"] |
| 11 | + data = delete!(pm, "user_defined_params") |
| 12 | + else |
| 13 | + data = pm |
| 14 | + end |
| 15 | + |
| 16 | + # additional modification |
| 17 | + for (i, branch) in pm["branch"] |
| 18 | + pm["branch"][i]["angmin"] = -6.28 |
| 19 | + pm["branch"][i]["angmax"] = 6.28 |
| 20 | + end |
| 21 | + |
| 22 | + # use build_ref to filter out inactive components |
| 23 | + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] |
| 24 | + |
| 25 | + ############################################################################### |
| 26 | + # 1. Building the Optimal Power Flow Model |
| 27 | + ############################################################################### |
| 28 | + # Initialize a JuMP Optimization Model |
| 29 | + #------------------------------------- |
| 30 | + model = JuMP.Model(solver) |
| 31 | + |
| 32 | + # Add voltage angles va for each bus |
| 33 | + JuMP.@variable(model, va[i in keys(ref[:bus])]) |
| 34 | + # note: [i in keys(ref[:bus])] adds one `va` variable for each bus in the network |
| 35 | + |
| 36 | + # Add voltage angles vm for each bus |
| 37 | + JuMP.@variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0) |
| 38 | + # note: this vairable also includes the voltage magnitude limits and a starting value |
| 39 | + |
| 40 | + # Add active power generation variable pg for each generator (including limits) |
| 41 | + JuMP.@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"], start=ref[:gen][i]["pg"]) |
| 42 | + # Add reactive power generation variable qg for each generator (including limits) |
| 43 | + JuMP.@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"], start=ref[:gen][i]["qg"]) |
| 44 | + |
| 45 | + # Add power flow variables p to represent the active power flow for each branch |
| 46 | + JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]) |
| 47 | + # Add power flow variables q to represent the reactive power flow for each branch |
| 48 | + JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]) |
| 49 | + # note: ref[:arcs] includes both the from (i,j) and the to (j,i) sides of a branch |
| 50 | + |
| 51 | + # Add Objective Function |
| 52 | + # ---------------------- |
| 53 | + |
| 54 | + JuMP.@objective(model, Min, sum((q[(content["element_index"], |
| 55 | + content["f_bus"], |
| 56 | + content["t_bus"])] - content["value"])^2 |
| 57 | + for (i, content) in user_defined_params["setpoint_q"])) |
| 58 | + |
| 59 | + # Add Constraints |
| 60 | + for (i,bus) in ref[:ref_buses] |
| 61 | + JuMP.@constraint(model, va[i] == 0) |
| 62 | + end |
| 63 | + |
| 64 | + # Nodal power balance constraints |
| 65 | + for (i,bus) in ref[:bus] |
| 66 | + # Build a list of the loads and shunt elements connected to the bus i |
| 67 | + bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]] |
| 68 | + bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]] |
| 69 | + |
| 70 | + # Active power balance at node i |
| 71 | + JuMP.@constraint(model, |
| 72 | + sum(p[a] for a in ref[:bus_arcs][i]) + # sum of active power flow on lines from bus i + |
| 73 | + sum(p_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) == # sum of active power flow on HVDC lines from bus i = |
| 74 | + sum(pg[g] for g in ref[:bus_gens][i]) - # sum of active power generation at bus i - |
| 75 | + sum(load["pd"] for load in bus_loads) - # sum of active load consumption at bus i - |
| 76 | + sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2 # sum of active shunt element injections at bus i |
| 77 | + ) |
| 78 | + |
| 79 | + # Reactive power balance at node i |
| 80 | + JuMP.@constraint(model, |
| 81 | + sum(q[a] for a in ref[:bus_arcs][i]) + # sum of reactive power flow on lines from bus i + |
| 82 | + sum(q_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) == # sum of reactive power flow on HVDC lines from bus i = |
| 83 | + sum(qg[g] for g in ref[:bus_gens][i]) - # sum of reactive power generation at bus i - |
| 84 | + sum(load["qd"] for load in bus_loads) + # sum of reactive load consumption at bus i - |
| 85 | + sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2 # sum of reactive shunt element injections at bus i |
| 86 | + ) |
| 87 | + end |
| 88 | + |
| 89 | + # Branch power flow physics and limit constraints |
| 90 | + for (i,branch) in ref[:branch] |
| 91 | + # Build the from variable id of the i-th branch, which is a tuple given by (branch id, from bus, to bus) |
| 92 | + f_idx = (i, branch["f_bus"], branch["t_bus"]) |
| 93 | + # Build the to variable id of the i-th branch, which is a tuple given by (branch id, to bus, from bus) |
| 94 | + t_idx = (i, branch["t_bus"], branch["f_bus"]) |
| 95 | + # note: it is necessary to distinguish between the from and to sides of a branch due to power losses |
| 96 | + |
| 97 | + p_fr = p[f_idx] # p_fr is a reference to the optimization variable p[f_idx] |
| 98 | + q_fr = q[f_idx] # q_fr is a reference to the optimization variable q[f_idx] |
| 99 | + p_to = p[t_idx] # p_to is a reference to the optimization variable p[t_idx] |
| 100 | + q_to = q[t_idx] # q_to is a reference to the optimization variable q[t_idx] |
| 101 | + # note: adding constraints to p_fr is equivalent to adding constraints to p[f_idx], and so on |
| 102 | + |
| 103 | + vm_fr = vm[branch["f_bus"]] # vm_fr is a reference to the optimization variable vm on the from side of the branch |
| 104 | + vm_to = vm[branch["t_bus"]] # vm_to is a reference to the optimization variable vm on the to side of the branch |
| 105 | + va_fr = va[branch["f_bus"]] # va_fr is a reference to the optimization variable va on the from side of the branch |
| 106 | + va_to = va[branch["t_bus"]] # va_fr is a reference to the optimization variable va on the to side of the branch |
| 107 | + |
| 108 | + # Compute the branch parameters and transformer ratios from the data |
| 109 | + g, b = PowerModels.calc_branch_y(branch) |
| 110 | + tr, ti = PowerModels.calc_branch_t(branch) |
| 111 | + g_fr = branch["g_fr"] |
| 112 | + b_fr = branch["b_fr"] |
| 113 | + g_to = branch["g_to"] |
| 114 | + b_to = branch["b_to"] |
| 115 | + tm = branch["tap"]^2 |
| 116 | + # note: tap is assumed to be 1.0 on non-transformer branches |
| 117 | + |
| 118 | + |
| 119 | + # AC Power Flow Constraints |
| 120 | + |
| 121 | + # From side of the branch flow |
| 122 | + JuMP.@NLconstraint(model, p_fr == (g+g_fr)/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) ) |
| 123 | + JuMP.@NLconstraint(model, q_fr == -(b+b_fr)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) ) |
| 124 | + |
| 125 | + # To side of the branch flow |
| 126 | + JuMP.@NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) ) |
| 127 | + JuMP.@NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) ) |
| 128 | + |
| 129 | + # Voltage angle difference limit |
| 130 | + JuMP.@constraint(model, va_fr - va_to <= branch["angmax"]) |
| 131 | + JuMP.@constraint(model, va_fr - va_to >= branch["angmin"]) |
| 132 | + |
| 133 | + # Apparent power limit, from side and to side |
| 134 | + JuMP.@constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2) |
| 135 | + JuMP.@constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2) |
| 136 | + end |
| 137 | + |
| 138 | + |
| 139 | + ############################################################################### |
| 140 | + # 3. Solve the Optimal Power Flow Model and Review the Results |
| 141 | + ############################################################################### |
| 142 | + JuMP.optimize!(model) |
| 143 | + |
| 144 | + ############################################################################### |
| 145 | + # 4. Create Result Dictionary such that the PowerModels Results can be used by pandapower |
| 146 | + ############################################################################### |
| 147 | + solution = Dict{String,Any}() |
| 148 | + push!(solution, "baseMVA" => data["baseMVA"]) |
| 149 | + push!(solution, "per_unit" => data["per_unit"]) |
| 150 | + push!(solution, "gen" => data["gen"]) |
| 151 | + push!(solution, "bus" => data["bus"]) |
| 152 | + push!(solution, "branch" => data["branch"]) |
| 153 | + push!(solution, "multinetwork" => false) |
| 154 | + push!(solution, "multiinfrastructrue" => false) |
| 155 | + |
| 156 | + for (i, gen) in solution["gen"] |
| 157 | + index = gen["index"] |
| 158 | + gen["qg"] = value(model[:qg][index]) |
| 159 | + gen["pg"] = value(model[:pg][index]) |
| 160 | + end |
| 161 | + |
| 162 | + for (i, bus) in solution["bus"] |
| 163 | + index = bus["index"] |
| 164 | + bus["vm"] = value(model[:vm][index]) |
| 165 | + bus["va"] = value(model[:va][index]) |
| 166 | + end |
| 167 | + |
| 168 | + for (i, branch) in solution["branch"] |
| 169 | + index = branch["index"] |
| 170 | + push!(branch, "qf" => value(model[:q][(index, branch["f_bus"], branch["t_bus"])])) |
| 171 | + push!(branch, "qt" => value(model[:q][(index, branch["t_bus"], branch["f_bus"])])) |
| 172 | + push!(branch, "pf" => value(model[:p][(index, branch["f_bus"], branch["t_bus"])])) |
| 173 | + push!(branch, "pt" => value(model[:p][(index, branch["t_bus"], branch["f_bus"])])) |
| 174 | + end |
| 175 | + println("The solver termination status is") |
| 176 | + |
| 177 | + result = Dict{String,Any}( |
| 178 | + "optimizer" => JuMP.solver_name(model), |
| 179 | + "termination_status" => JuMP.termination_status(model), # "LOCALLY_SOLVED", |
| 180 | + "primal_status" => JuMP.primal_status(model), |
| 181 | + "dual_status" => JuMP.dual_status(model), |
| 182 | + "objective" => InfrastructureModels._guard_objective_value(model), |
| 183 | + "objective_lb" => InfrastructureModels._guard_objective_bound(model), |
| 184 | + "solve_time" => solve_time, |
| 185 | + "solution" => solution) |
| 186 | + |
| 187 | + return result |
| 188 | +end |
| 189 | + |
| 190 | + |
| 191 | +# json_path = "C:/Users/fmeier/pandapower/pandapower/test/opf/case5_clm_matfile_va.json" |
| 192 | +# # #@enter run_powermodels(json_path) |
| 193 | +# # |
| 194 | +# result = run_powermodels(json_path) |
| 195 | +# println(result["termination_status"] == LOCALLY_SOLVED) |
| 196 | +# println(isapprox(result["objective"], 17015.5; atol = 1e0)) |
| 197 | +# mit eingeschränkter slack spannung: 17082.819507648066 |
0 commit comments