diff --git a/src/GridCalEngine/Simulations/OPF/NumericalMethods/ac_opf.py b/src/GridCalEngine/Simulations/OPF/NumericalMethods/ac_opf.py index 7662e5db6..1b6231070 100644 --- a/src/GridCalEngine/Simulations/OPF/NumericalMethods/ac_opf.py +++ b/src/GridCalEngine/Simulations/OPF/NumericalMethods/ac_opf.py @@ -620,9 +620,8 @@ def ac_optimal_power_flow(nc: NumericalCircuit, if opf_options.acopf_mode == AcOpfMode.ACOPFslacks: nsl = 2 * npq + 2 * n_br_mon # Slack relaxations for constraints - c_s = 1 * np.power(nc.branch_data.overload_cost[br_mon_idx] + 0.1, - 1.0) # Cost squared since the slack is also squared - c_v = 1 * (nc.bus_data.cost_v[pq] + 0.1) + c_s = 1 * np.power(nc.branch_data.overload_cost[br_mon_idx] + 0.1, 1.0) # Cost squared since the slack is also squared + c_v = 1000 * (nc.bus_data.cost_v[pq] + 0.1) sl_sf0 = np.ones(n_br_mon) sl_st0 = np.ones(n_br_mon) sl_vmax0 = np.ones(npq) @@ -778,6 +777,14 @@ def ac_optimal_power_flow(nc: NumericalCircuit, Sf = result.structs.Sf St = result.structs.St loading = np.abs(Sf) / (rates + 1e-9) + + if opf_options.acopf_mode == AcOpfMode.ACOPFslacks: + overloads_sf = (np.power(np.power(rates, 2) + sl_sf, 0.5) - rates)*Sbase + overloads_st = (np.power(np.power(rates, 2) + sl_st, 0.5) - rates)*Sbase + + else: + pass + hvdc_power = nc.hvdc_data.Pset.copy() hvdc_power[hvdc_disp_idx] = Pfdc hvdc_loading = hvdc_power / (nc.hvdc_data.rate + 1e-9) @@ -787,9 +794,7 @@ def ac_optimal_power_flow(nc: NumericalCircuit, tap_phase[k_tau] = tapt Pcost = np.zeros(nc.ngen) Pcost[gen_disp_idx] = c0 + c1 * Pg[gen_disp_idx] + c2 * np.power(Pg[gen_disp_idx], 2.0) - Pcost[gen_nondisp_idx] = c0n + c1n * np.real(Sg_undis) + c2n * np.power(np.real(Sg_undis), 2.0) - nodal_capacity = slcap * Sbase if opf_options.verbose > 0: @@ -885,19 +890,19 @@ def ac_optimal_power_flow(nc: NumericalCircuit, if opf_options.acopf_mode == AcOpfMode.ACOPFslacks: for k in range(n_br_mon): - if sl_sf[k] > opf_options.ips_tolerance: - logger.add_warning('Branch overload in the from sense', + if overloads_sf[k] > opf_options.ips_tolerance * Sbase: + logger.add_warning('Branch overload in the from sense (MVA)', device=str(br_mon_idx[k]), device_property="Slack", - value=str(sl_sf[k]), - expected_value=f'< {opf_options.ips_tolerance}') + value=str(overloads_sf[k]), + expected_value=f'< {opf_options.ips_tolerance*Sbase}') - if sl_st[k] > opf_options.ips_tolerance: - logger.add_warning('Branch overload in the to sense', + if overloads_st[k] > opf_options.ips_tolerance * Sbase: + logger.add_warning('Branch overload in the to sense (MVA)', device=str(br_mon_idx[k]), device_property="Slack", - value=str(sl_st[k]), - expected_value=f'< {opf_options.ips_tolerance}') + value=str(overloads_st[k]), + expected_value=f'< {opf_options.ips_tolerance*Sbase}') for i in range(npq): if sl_vmax[i] > opf_options.ips_tolerance: diff --git a/src/GridCalEngine/Utils/NumericalMethods/ips.py b/src/GridCalEngine/Utils/NumericalMethods/ips.py index b4bdf5cb4..8ce6ca9c6 100644 --- a/src/GridCalEngine/Utils/NumericalMethods/ips.py +++ b/src/GridCalEngine/Utils/NumericalMethods/ips.py @@ -334,7 +334,7 @@ def interior_point_solver(x0: Vec, feascond = calc_feascond(g=ret.G, h=ret.H, x=x, z=z) converged = error <= gamma - + maxdispl = 0 error_evolution = np.zeros(max_iter + 1) feascond_evolution = np.zeros(max_iter + 1) @@ -414,7 +414,7 @@ def interior_point_solver(x0: Vec, feascond = max([g_norm, max(ret.H)]) / (1 + max([np.linalg.norm(x, np.Inf), z_norm])) gradcond = np.linalg.norm(lx, np.Inf) / (1 + max([lam_norm, mu_norm])) error = np.max([feascond, gradcond, gamma]) - + maxdispl = np.max(np.r_[dx, dlam, dz, dmu]) z_inv = diags(1.0 / z) mu_diag = diags(mu) @@ -432,6 +432,7 @@ def interior_point_solver(x0: Vec, print("INEQ:\n", ineq_df) print("\tGamma:", gamma) print("\tErr:", error) + print("\tMax Displacement:", maxdispl) # Add an iteration step iter_counter += 1 @@ -449,6 +450,7 @@ def interior_point_solver(x0: Vec, print(f"\tF.obj: {ret.f * 1e4}") print(f"\tErr: {error}") print(f'\tIterations: {iter_counter}') + print(f'\tMax Displacement: {maxdispl}' ) print(f'\tTime elapsed (s): {t_end - t_start}') print(f'\tFeas cond: ', feascond) diff --git a/src/trunk/acopf/acopf_run.py b/src/trunk/acopf/acopf_run.py index 791862a5c..7efc38b1b 100644 --- a/src/trunk/acopf/acopf_run.py +++ b/src/trunk/acopf/acopf_run.py @@ -8,6 +8,7 @@ from GridCalEngine.Simulations.NodalCapacity.nodal_capacity_options import NodalCapacityOptions import numpy as np import pandas as pd +import math from GridCalEngine.enumerations import NodalCapacityMethod @@ -398,21 +399,31 @@ def casehvdc(): cwd = os.getcwd() # Go back two directories - new_directory = os.path.abspath(os.path.join(cwd, '..', '..', '..')) + new_directory = os.path.abspath(os.path.join(cwd, '..', '..', '..', '..')) - file_path = os.path.join(new_directory, 'Grids_and_profiles', 'grids', 'entrada_a_aopf.raw') + file_path = os.path.join(new_directory, 'REE Grids', 'entrada_a_aopf.raw') grid = gce.FileOpen(file_path).open() + for gen in grid.generators: + gen.qmin_set = -0.8 * gen.Snom + gen.qmax_set = 0.8 * gen.Snom + options = gce.PowerFlowOptions(gce.SolverType.NR, verbose=False) power_flow = gce.PowerFlowDriver(grid, options) power_flow.run() - opf_options = gce.OptimalPowerFlowOptions(solver=gce.SolverType.NONLINEAR_OPF, acopf_mode=AcOpfMode.ACOPFslacks, + opf_options = gce.OptimalPowerFlowOptions(solver=gce.SolverType.NONLINEAR_OPF, acopf_mode=AcOpfMode.ACOPFstd, verbose=1, ips_iterations=150, ips_tolerance=1e-8) pf_options = gce.PowerFlowOptions(solver_type=gce.SolverType.NR, verbose=3) - run_nonlinear_opf(grid=grid, pf_options=pf_options, opf_options=opf_options, plot_error=True, pf_init=True) - + opf_options = gce.OptimalPowerFlowOptions(solver=gce.SolverType.NONLINEAR_OPF, ips_tolerance=1e-8, + ips_iterations=50, verbose=1, acopf_mode=AcOpfMode.ACOPFslacks) + res = run_nonlinear_opf(grid=grid, pf_options=pf_options, opf_options=opf_options, plot_error=True, pf_init=True, + optimize_nodal_capacity=False, + nodal_capacity_sign=-1.0, + capacity_nodes_idx=np.array([10])) + #run_nonlinear_opf(grid=grid, pf_options=pf_options, opf_options=opf_options, plot_error=True, pf_init=True) + print('') def caseREE(): """ @@ -423,11 +434,15 @@ def caseREE(): # Go back two directories new_directory = os.path.abspath(os.path.join(cwd, '..', '..', '..', '..')) - # file_path = os.path.join(new_directory, 'REE Grids', 'entrada_a_aopf.raw') - file_path = 'C:/Users/J/Documents/ree_opf/entrada_a_aopf.raw' + file_path = os.path.join(new_directory, 'REE Grids', 'entrada_a_aopf.raw') + # file_path = 'C:/Users/J/Documents/ree_opf/entrada_a_aopf.raw' grid = gce.FileOpen(file_path).open() + for gen in grid.generators: + gen.qmin_set = -0.8 * gen.Snom + gen.qmax_set = 0.8 * gen.Snom + disp_areas = ['A11', 'A15'] dict_bus_lims = {'21215': [230, 225], '11055': [410, 405], @@ -436,20 +451,23 @@ def caseREE(): '15005': [410, 405], '15015': [410, 405]} tol = 1e-4 - vm_cost = 1e2 + vm_cost = 1e4 i = 0 for gen in grid.generators: if gen.bus.area.name in disp_areas: # P limits -> restrict them very close to P print(f'Select generator {i}') - gen.Pmax = gen.P + tol - gen.Pmin = gen.P - tol + #gen.Pmax = gen.P #+ tol + #gen.Pmin = gen.P #- tol # Tanmax -> set pf close to 0 to get large tanmax - gen.Pf = tol + #gen.Pf = tol else: - gen.enabled_dispatch = False - gen.Pmax = gen.P + tol - gen.Pmin = gen.P - tol + + #gen.enabled_dispatch = False + gen.Pmax = gen.P #+ tol + gen.Pmin = gen.P #- tol + #gen.Qmax = abs(math.tan(math.acos(gen.Pf)) * gen.P) + #gen.Qmin = -abs(math.tan(math.acos(gen.Pf)) * gen.P) # i += 1 @@ -472,7 +490,7 @@ def caseREE(): power_flow = gce.PowerFlowDriver(grid, options) power_flow.run() - opf_options = gce.OptimalPowerFlowOptions(solver=gce.SolverType.NONLINEAR_OPF, acopf_mode=AcOpfMode.ACOPFstd, + opf_options = gce.OptimalPowerFlowOptions(solver=gce.SolverType.NONLINEAR_OPF, acopf_mode=AcOpfMode.ACOPFslacks, verbose=1, ips_iterations=100, ips_tolerance=1e-8) pf_options = gce.PowerFlowOptions(solver_type=gce.SolverType.NR, verbose=3) run_nonlinear_opf(grid=grid, pf_options=pf_options, opf_options=opf_options, plot_error=True, pf_init=True) @@ -508,13 +526,13 @@ def case_nodalcap(): # linn5bus_example() # two_grids_of_3bus() # case9() - case14_linear_vs_nonlinear() + #case14_linear_vs_nonlinear() # case14() # case_gb() # case6ww() # case_pegase89() # case300() # casepegase13k() - # casehvdc() - # caseREE() + #casehvdc() + caseREE() # case_nodalcap()