From 0cc3757416a4fae8a628eef7151e62db8680ea8d Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Fri, 14 May 2021 16:06:26 -0600 Subject: [PATCH 1/6] pw linear cost support in acopf --- egret/model_library/transmission/gen.py | 88 +++++++++++++++++++++---- egret/models/acopf.py | 23 +++++-- 2 files changed, 93 insertions(+), 18 deletions(-) diff --git a/egret/model_library/transmission/gen.py b/egret/model_library/transmission/gen.py index fbd0ce05..34be3abb 100644 --- a/egret/model_library/transmission/gen.py +++ b/egret/model_library/transmission/gen.py @@ -28,6 +28,65 @@ def declare_var_qg(model, index_set, **kwargs): decl.declare_var('qg', model=model, index_set=index_set, **kwargs) +def declare_var_pg_cost(model, index_set, **kwargs): + decl.declare_var('pg_cost', model, index_set=index_set, **kwargs) + + +def declare_var_qg_cost(model, index_set, **kwargs): + decl.declare_var('qg_cost', model, index_set=index_set, **kwargs) + + +def _pw_cost_helper(cost_dict, cost_var, gen_var, pw_cost_set, gen_name, indexed_pw_cost_con): + if cost_dict['cost_curve_type'] == 'polynomial': + pass + elif cost_dict['cost_curve_type'] == 'piecewise': + pt0, cost0 = cost_dict['values'][0] + if gen_var.lb < pt0: + raise ValueError('Piecewise costs require that the lower bound on pg be greater than or equal to the first piecewise point.') + + last_slope = None + for ndx in range(len(cost_dict['values']) - 1): + pt1, cost1 = cost_dict['values'][ndx] + pt2, cost2 = cost_dict['values'][ndx + 1] + slope = (cost2 - cost1) / (pt2 - pt1) + if last_slope is not None: + if slope < last_slope: + all_slopes = [(cost_dict['values'][i + 1][1] - cost_dict['values'][i][1]) / (cost_dict['values'][i + 1][0] - cost_dict['values'][i][0]) for i in range(len(cost_dict['values']) - 1)] + raise ValueError(f'Piecewise costs must be convex; generator: {gen_name}; slopes: {all_slopes}') + intercept = cost2 - slope * pt2 + pw_cost_set.add((gen_name, ndx)) + indexed_pw_cost_con[gen_name, ndx] = cost_var >= slope * gen_var + intercept + last_slope = slope + else: + raise ValueError(f"Unrecognized cost_cureve_type: {cost_dict['cost_curve_type']}") + + +def declare_piecewise_cost_cons(model, index_set, p_costs, q_costs=None): + m = model + + m.pg_piecewise_cost_set = pe.Set(dimen=2) + m.qg_piecewise_cost_set = pe.Set(dimen=2) + m.pg_piecewise_cost_cons = pe.Constraint(m.pg_piecewise_cost_set) + m.qg_piecewise_cost_cons = pe.Constraint(m.qg_piecewise_cost_set) + + for gen_name in index_set: + if p_costs is not None and gen_name in p_costs: + _pw_cost_helper(cost_dict=p_costs[gen_name], + cost_var=m.pg_cost[gen_name], + gen_var=m.pg[gen_name], + pw_cost_set=m.pg_piecewise_cost_set, + gen_name=gen_name, + indexed_pw_cost_con=m.pg_piecewise_cost_cons) + + if q_costs is not None and gen_name in q_costs: + _pw_cost_helper(cost_dict=q_costs[gen_name], + cost_var=m.qg_cost[gen_name], + gen_var=m.qg[gen_name], + pw_cost_set=m.qg_piecewise_cost_set, + gen_name=gen_name, + indexed_pw_cost_con=m.qg_piecewise_cost_cons) + + def declare_expression_pgqg_operating_cost(model, index_set, p_costs, q_costs=None): """ @@ -41,20 +100,23 @@ def declare_expression_pgqg_operating_cost(model, index_set, m.pg_operating_cost = pe.Expression(expr_set) m.qg_operating_cost = pe.Expression(expr_set) - found_q_costs = False for gen_name in expr_set: if p_costs is not None and gen_name in p_costs: - #TODO: We assume that the costs are polynomial type - assert p_costs[gen_name]['cost_curve_type'] == 'polynomial' - m.pg_operating_cost[gen_name] = \ - sum(v*m.pg[gen_name]**i for i, v in p_costs[gen_name]['values'].items()) + if p_costs[gen_name]['cost_curve_type'] == 'polynomial': + m.pg_operating_cost[gen_name] = sum(v*m.pg[gen_name]**i for i, v in p_costs[gen_name]['values'].items()) + elif p_costs[gen_name]['cost_curve_type'] == 'piecewise': + m.pg_operating_cost[gen_name] = m.pg_cost[gen_name] + else: + raise ValueError(f"Unrecognized cost_cureve_type: {p_costs[gen_name]['cost_curve_type']}") + else: + m.pg_operating_cost[gen_name] = 0 if q_costs is not None and gen_name in q_costs: - #TODO: We assume that the costs are polynomial type - assert q_costs['cost_curve_type'] == 'polynomial' - found_q_costs = True - m.qg_operating_cost[gen_name] = \ - sum(v*m.qg[gen_name]**i for i, v in q_costs[gen_name]['values'].items()) - - if not found_q_costs: - del m.qg_operating_cost + if q_costs[gen_name]['cost_curve_type'] == 'polynomial': + m.qg_operating_cost[gen_name] = sum(v*m.qg[gen_name]**i for i, v in q_costs[gen_name]['values'].items()) + elif q_costs[gen_name]['cost_curve_type'] == 'piecewise': + m.qg_operating_cost[gen_name] = m.qg_cost[gen_name] + else: + raise ValueError(f"Unrecognized cost_cureve_type: {q_costs[gen_name]['cost_curve_type']}") + else: + m.qg_operating_cost[gen_name] = 0 diff --git a/egret/models/acopf.py b/egret/models/acopf.py index ca16f464..530fc69a 100644 --- a/egret/models/acopf.py +++ b/egret/models/acopf.py @@ -228,18 +228,31 @@ def _create_base_power_ac_model(model_data, include_feasibility_slack=False): branches=branches ) - ### declare the generator cost objective + # declare the generator cost objective + pw_cost_gens = list() + p_cost = gen_attrs['p_cost'] + if 'q_cost' in gen_attrs: + q_cost = gen_attrs['q_cost'] + else: + q_cost = dict() + for gen_name in gen_attrs['names']: + if p_cost[gen_name]['cost_curve_type'] == 'piecewise' or (gen_name in q_cost and q_cost[gen_name]['cost_curve_type'] == 'piecewise'): + pw_cost_gens.append(gen_name) + + if len(pw_cost_gens) > 0: + libgen.declare_var_pg_cost(model=model, index_set=pw_cost_gens) + libgen.declare_var_qg_cost(model=model, index_set=pw_cost_gens) + libgen.declare_piecewise_cost_cons(model=model, index_set=pw_cost_gens, p_costs=p_cost, q_costs=q_cost) + libgen.declare_expression_pgqg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=gen_attrs['p_cost'], - q_costs=gen_attrs.get('q_cost', None) - ) + q_costs=gen_attrs.get('q_cost', None)) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + obj_expr += sum(model.qg_operating_cost[gen_name] for gen_name in model.qg_operating_cost) if include_feasibility_slack: obj_expr += penalty_expr - if hasattr(model, 'qg_operating_cost'): - obj_expr += sum(model.qg_operating_cost[gen_name] for gen_name in model.qg_operating_cost) model.obj = pe.Objective(expr=obj_expr) From 738e4e1199e6220751b75b481459518a9546ed47 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Sat, 15 May 2021 03:15:33 -0600 Subject: [PATCH 2/6] pw linear cost support in acopf --- egret/models/tests/test_acopf.py | 49 ++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/egret/models/tests/test_acopf.py b/egret/models/tests/test_acopf.py index 64157243..b189bf0d 100644 --- a/egret/models/tests/test_acopf.py +++ b/egret/models/tests/test_acopf.py @@ -20,6 +20,7 @@ from egret.parsers.matpower_parser import create_ModelData import numpy.testing as npt import json +import numpy as np current_dir = os.path.dirname(os.path.abspath(__file__)) case_names = ['pglib_opf_case3_lmbd','pglib_opf_case30_ieee','pglib_opf_case300_ieee','pglib_opf_case3012wp_k'] @@ -162,5 +163,53 @@ def test_acopf_model(self, test_case, soln_case, p_and_v_soln_case, include_kwar _test_p_and_v(self, p_and_v_soln_case, md) +def poly_cost_to_pw_cost(md: ModelData, num_points=10): + gen_attrs = md.attributes(element_type='generator') + p_cost = gen_attrs['p_cost'] + for gen_name, cost_dict in p_cost.items(): + assert cost_dict['cost_curve_type'] == 'polynomial' + cost_dict['cost_curve_type'] = 'piecewise' + poly_coefs = cost_dict['values'] + pw_values = list() + p_min = gen_attrs['p_min'][gen_name] + p_max = gen_attrs['p_max'][gen_name] + for pt in np.linspace(p_min, p_max, num_points): + cost = sum(coef*pt**exponent for exponent, coef in poly_coefs.items()) + pw_values.append((pt, cost)) + cost_dict['values'] = pw_values + + +class TestPWCost(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + download_dir = os.path.join(current_dir, 'transmission_test_instances') + if not os.path.exists(os.path.join(download_dir, 'pglib-opf-master')): + from egret.thirdparty.get_pglib_opf import get_pglib_opf + get_pglib_opf(download_dir) + + def test_case30_pw_cost(self): + original_md = ModelData.read(os.path.join(current_dir, 'transmission_test_instances', 'pglib-opf-master', 'pglib_opf_case30_as.m')) + + md = original_md.clone() + m_poly, scaled_md = create_atan_acopf_model(md) + opt = pe.SolverFactory('ipopt') + res = opt.solve(m_poly) + poly_obj = pe.value(m_poly.obj.expr) + self.assertAlmostEqual(poly_obj, 803.12692652061719, places=2) + + poly_cost_to_pw_cost(md, num_points=2) + m_pw, scaled_md = create_atan_acopf_model(md) + res = opt.solve(m_pw) + pw_obj = pe.value(m_pw.obj.expr) + self.assertAlmostEqual(pw_obj, 827.62708294193681, places=2) + + md = original_md.clone() + poly_cost_to_pw_cost(md, num_points=10) + m_pw, scaled_md = create_atan_acopf_model(md) + res = opt.solve(m_pw) + pw_obj = pe.value(m_pw.obj.expr) + self.assertAlmostEqual(pw_obj, 803.56080829371604, places=2) + + if __name__ == '__main__': unittest.main() From 8af042ae80f5778836a21a979a67b855acb71fe0 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 17 May 2021 15:20:57 -0600 Subject: [PATCH 3/6] support for pw costs in acopf and dcopf --- egret/model_library/transmission/gen.py | 238 ++++++++++++++---- .../transmission/tests/test_tx_utils.py | 202 +++++++++++++++ egret/model_library/transmission/tx_utils.py | 103 ++++++++ egret/model_library/unit_commitment/params.py | 71 ++---- egret/models/acopf.py | 46 ++-- egret/models/dcopf.py | 40 +-- egret/models/tests/test_dcopf.py | 82 ++++++ 7 files changed, 642 insertions(+), 140 deletions(-) create mode 100644 egret/model_library/transmission/tests/test_tx_utils.py diff --git a/egret/model_library/transmission/gen.py b/egret/model_library/transmission/gen.py index 34be3abb..f0659fe7 100644 --- a/egret/model_library/transmission/gen.py +++ b/egret/model_library/transmission/gen.py @@ -13,6 +13,9 @@ """ import pyomo.environ as pe import egret.model_library.decl as decl +from egret.model_library.transmission import tx_utils +from pyomo.core.expr.numeric_expr import LinearExpression + def declare_var_pg(model, index_set, **kwargs): """ @@ -28,94 +31,235 @@ def declare_var_qg(model, index_set, **kwargs): decl.declare_var('qg', model=model, index_set=index_set, **kwargs) -def declare_var_pg_cost(model, index_set, **kwargs): - decl.declare_var('pg_cost', model, index_set=index_set, **kwargs) +def pw_gen_generator(index_set, costs): + for gen_name in index_set: + if gen_name not in costs: + continue + curve = costs[gen_name] + assert curve['cost_curve_type'] in {'piecewise', 'polynomial'} + if curve['cost_curve_type'] != 'piecewise': + continue + yield gen_name + + +def declare_var_delta_pg(model, index_set, p_costs): + m = model + + m.delta_pg_set = pe.Set(dimen=2) + m.delta_pg = pe.Var(m.delta_pg_set) + for gen_name in pw_gen_generator(index_set, p_costs): + p_min = m.pg[gen_name].lb + p_max = m.pg[gen_name].ub + curve = p_costs[gen_name] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=p_min, + p_max=p_max, + gen_name=gen_name) + for ndx, ((o1, c1), (o2, c2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): + m.delta_pg_set.add((gen_name, ndx)) + m.delta_pg[gen_name, ndx].setlb(0) + m.delta_pg[gen_name, ndx].setub(o2 - o1) + + +def declare_var_delta_qg(model, index_set, q_costs): + m = model + + m.delta_qg_set = pe.Set(dimen=2) + m.delta_qg = pe.Var(m.delta_qg_set) + for gen_name in pw_gen_generator(index_set, q_costs): + q_min = m.qg[gen_name].lb + q_max = m.qg[gen_name].ub + curve = q_costs[gen_name] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=q_min, + p_max=q_max, + gen_name=gen_name) + for ndx, ((o1, c1), (o2, c2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): + m.delta_qg_set.add((gen_name, ndx)) + m.delta_qg[gen_name, ndx].setlb(0) + m.delta_qg[gen_name, ndx].setub(o2 - o1) + + +def declare_pg_delta_pg_con(model, index_set, p_costs): + m = model + + m.pg_delta_pg_con_set = pe.Set() + m.pg_delta_pg_con = pe.Constraint(m.pg_delta_pg_con_set) + for gen_name in pw_gen_generator(index_set, p_costs): + p_min = m.pg[gen_name].lb + p_max = m.pg[gen_name].ub + curve = p_costs[gen_name] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=p_min, + p_max=p_max, + gen_name=gen_name) + m.pg_delta_pg_con_set.add(gen_name) + lin_coefs = [-1] + lin_vars = [m.pg[gen_name]] + for ndx, ((o1, c1), (o2, c2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): + lin_coefs.append(1) + lin_vars.append(m.delta_pg[gen_name, ndx]) + expr = LinearExpression(constant=cleaned_values[0][0], linear_coefs=lin_coefs, linear_vars=lin_vars) + m.pg_delta_pg_con[gen_name] = (expr, 0) + + +def declare_qg_delta_qg_con(model, index_set, q_costs): + m = model + + m.qg_delta_qg_con_set = pe.Set() + m.qg_delta_qg_con = pe.Constraint(m.qg_delta_qg_con_set) + for gen_name in pw_gen_generator(index_set, q_costs): + q_min = m.qg[gen_name].lb + q_max = m.qg[gen_name].ub + curve = q_costs[gen_name] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=q_min, + p_max=q_max, + gen_name=gen_name) + m.qg_delta_qg_con_set.add(gen_name) + lin_coefs = [-1] + lin_vars = [m.qg[gen_name]] + for ndx, ((o1, c1), (o2, c2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): + lin_coefs.append(1) + lin_vars.append(m.delta_qg[gen_name, ndx]) + expr = LinearExpression(constant=cleaned_values[0][0], linear_coefs=lin_coefs, linear_vars=lin_vars) + m.qg_delta_qg_con[gen_name] = (expr, 0) + + +def declare_var_pg_cost(model, index_set, p_costs): + m = model + actual_indices = list(pw_gen_generator(index_set, p_costs)) + m.pg_cost_set = pe.Set(initialize=actual_indices) + m.pg_cost = pe.Var(m.pg_cost_set) -def declare_var_qg_cost(model, index_set, **kwargs): - decl.declare_var('qg_cost', model, index_set=index_set, **kwargs) +def declare_var_qg_cost(model, index_set, q_costs): + m = model + actual_indices = list(pw_gen_generator(index_set, q_costs)) + m.qg_cost_set = pe.Set(initialize=actual_indices) + m.qg_cost = pe.Var(m.qg_cost_set) def _pw_cost_helper(cost_dict, cost_var, gen_var, pw_cost_set, gen_name, indexed_pw_cost_con): if cost_dict['cost_curve_type'] == 'polynomial': pass elif cost_dict['cost_curve_type'] == 'piecewise': - pt0, cost0 = cost_dict['values'][0] - if gen_var.lb < pt0: - raise ValueError('Piecewise costs require that the lower bound on pg be greater than or equal to the first piecewise point.') - - last_slope = None - for ndx in range(len(cost_dict['values']) - 1): - pt1, cost1 = cost_dict['values'][ndx] - pt2, cost2 = cost_dict['values'][ndx + 1] + cleaned_values = tx_utils.validate_and_clean_cost_curve(cost_dict, + curve_type='cost_curve', + p_min=gen_var.lb, + p_max=gen_var.ub, + gen_name=gen_name) + for ndx, ((pt1, cost1), (pt2, cost2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): slope = (cost2 - cost1) / (pt2 - pt1) - if last_slope is not None: - if slope < last_slope: - all_slopes = [(cost_dict['values'][i + 1][1] - cost_dict['values'][i][1]) / (cost_dict['values'][i + 1][0] - cost_dict['values'][i][0]) for i in range(len(cost_dict['values']) - 1)] - raise ValueError(f'Piecewise costs must be convex; generator: {gen_name}; slopes: {all_slopes}') intercept = cost2 - slope * pt2 pw_cost_set.add((gen_name, ndx)) indexed_pw_cost_con[gen_name, ndx] = cost_var >= slope * gen_var + intercept - last_slope = slope else: raise ValueError(f"Unrecognized cost_cureve_type: {cost_dict['cost_curve_type']}") -def declare_piecewise_cost_cons(model, index_set, p_costs, q_costs=None): +def declare_piecewise_pg_cost_cons(model, index_set, p_costs): m = model m.pg_piecewise_cost_set = pe.Set(dimen=2) - m.qg_piecewise_cost_set = pe.Set(dimen=2) m.pg_piecewise_cost_cons = pe.Constraint(m.pg_piecewise_cost_set) + + for gen_name in pw_gen_generator(index_set=index_set, costs=p_costs): + _pw_cost_helper(cost_dict=p_costs[gen_name], + cost_var=m.pg_cost[gen_name], + gen_var=m.pg[gen_name], + pw_cost_set=m.pg_piecewise_cost_set, + gen_name=gen_name, + indexed_pw_cost_con=m.pg_piecewise_cost_cons) + + +def declare_piecewise_qg_cost_cons(model, index_set, q_costs): + m = model + + m.qg_piecewise_cost_set = pe.Set(dimen=2) m.qg_piecewise_cost_cons = pe.Constraint(m.qg_piecewise_cost_set) - for gen_name in index_set: - if p_costs is not None and gen_name in p_costs: - _pw_cost_helper(cost_dict=p_costs[gen_name], - cost_var=m.pg_cost[gen_name], - gen_var=m.pg[gen_name], - pw_cost_set=m.pg_piecewise_cost_set, - gen_name=gen_name, - indexed_pw_cost_con=m.pg_piecewise_cost_cons) - - if q_costs is not None and gen_name in q_costs: - _pw_cost_helper(cost_dict=q_costs[gen_name], - cost_var=m.qg_cost[gen_name], - gen_var=m.qg[gen_name], - pw_cost_set=m.qg_piecewise_cost_set, - gen_name=gen_name, - indexed_pw_cost_con=m.qg_piecewise_cost_cons) - - -def declare_expression_pgqg_operating_cost(model, index_set, - p_costs, q_costs=None): + for gen_name in pw_gen_generator(index_set=index_set, costs=q_costs): + _pw_cost_helper(cost_dict=q_costs[gen_name], + cost_var=m.qg_cost[gen_name], + gen_var=m.qg[gen_name], + pw_cost_set=m.qg_piecewise_cost_set, + gen_name=gen_name, + indexed_pw_cost_con=m.qg_piecewise_cost_cons) + + +def declare_expression_pg_operating_cost(model, index_set, p_costs, pw_formulation='delta'): """ Create the Expression objects to represent the operating costs - for the real and reactive (if present) power of each of the - generators. + for the real power of each of the generators. """ m = model - expr_set = decl.declare_set('_expr_g_operating_cost', + expr_set = decl.declare_set('_expr_pg_operating_cost', model=model, index_set=index_set) m.pg_operating_cost = pe.Expression(expr_set) - m.qg_operating_cost = pe.Expression(expr_set) for gen_name in expr_set: - if p_costs is not None and gen_name in p_costs: + if gen_name in p_costs: if p_costs[gen_name]['cost_curve_type'] == 'polynomial': m.pg_operating_cost[gen_name] = sum(v*m.pg[gen_name]**i for i, v in p_costs[gen_name]['values'].items()) elif p_costs[gen_name]['cost_curve_type'] == 'piecewise': - m.pg_operating_cost[gen_name] = m.pg_cost[gen_name] + if pw_formulation == 'delta': + p_min = m.pg[gen_name].lb + p_max = m.pg[gen_name].ub + curve = p_costs[gen_name] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=p_min, + p_max=p_max, + gen_name=gen_name) + expr = cleaned_values[0][1] + for ndx, ((o1, c1), (o2, c2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): + slope = (c2 - c1) / (o2 - o1) + expr += slope * m.delta_pg[gen_name, ndx] + m.pg_operating_cost[gen_name] = expr + else: + m.pg_operating_cost[gen_name] = m.pg_cost[gen_name] else: raise ValueError(f"Unrecognized cost_cureve_type: {p_costs[gen_name]['cost_curve_type']}") else: m.pg_operating_cost[gen_name] = 0 - if q_costs is not None and gen_name in q_costs: + +def declare_expression_qg_operating_cost(model, index_set, q_costs, pw_formulation='delta'): + """ + Create the Expression objects to represent the operating costs + for the reactive power of each of the generators. + """ + m = model + expr_set = decl.declare_set('_expr_qg_operating_cost', + model=model, index_set=index_set) + m.qg_operating_cost = pe.Expression(expr_set) + + for gen_name in expr_set: + if gen_name in q_costs: if q_costs[gen_name]['cost_curve_type'] == 'polynomial': m.qg_operating_cost[gen_name] = sum(v*m.qg[gen_name]**i for i, v in q_costs[gen_name]['values'].items()) elif q_costs[gen_name]['cost_curve_type'] == 'piecewise': - m.qg_operating_cost[gen_name] = m.qg_cost[gen_name] + if pw_formulation == 'delta': + q_min = m.qg[gen_name].lb + q_max = m.qg[gen_name].ub + curve = q_costs[gen_name] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=q_min, + p_max=q_max, + gen_name=gen_name) + expr = cleaned_values[0][1] + for ndx, ((o1, c1), (o2, c2)) in enumerate(zip(cleaned_values, cleaned_values[1:])): + slope = (c2 - c1) / (o2 - o1) + expr += slope * m.delta_pg[gen_name, ndx] + m.qg_operating_cost[gen_name] = expr + else: + m.qg_operating_cost[gen_name] = m.qg_cost[gen_name] else: raise ValueError(f"Unrecognized cost_cureve_type: {q_costs[gen_name]['cost_curve_type']}") else: diff --git a/egret/model_library/transmission/tests/test_tx_utils.py b/egret/model_library/transmission/tests/test_tx_utils.py new file mode 100644 index 00000000..8e57bc5e --- /dev/null +++ b/egret/model_library/transmission/tests/test_tx_utils.py @@ -0,0 +1,202 @@ +import unittest +from egret.model_library.transmission import tx_utils +import logging +import copy + + +def _example_quadratic(p): + return 0.05 * p ** 2 + p + 3 + + +def example_pw_curve(): + curve = dict() + curve['data_type'] = 'cost_curve' + curve['cost_curve_type'] = 'piecewise' + curve['values'] = [(10, _example_quadratic(10)), + (30, _example_quadratic(30)), + (50, _example_quadratic(50)), + (70, _example_quadratic(70)), + (90, _example_quadratic(90))] + return curve + + +def example_poly_curve(): + curve = dict() + curve['data_type'] = 'cost_curve' + curve['cost_curve_type'] = 'polynomial' + curve['values'] = {0: 3, 1: 1, 2: 0.05} + return curve + + +class TestValidateCostCurves(unittest.TestCase): + def test_pw_simple(self): + curve = example_pw_curve() + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, curve['values']) + self.assertIsNot(cleaned_values, curve['values']) + + def test_wrong_curve_type(self): + curve = example_pw_curve() + curve['data_type'] = 'blah' + with self.assertRaises(ValueError): + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + + def test_pw_no_values(self): + curve = example_pw_curve() + curve['values'] = list() + with self.assertLogs('egret.model_library.transmission.tx_utils', level=logging.WARNING) as cm: + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, curve['values']) + self.assertIsNot(cleaned_values, curve['values']) + self.assertEqual(cm.output, ['WARNING:egret.model_library.transmission.tx_utils:WARNING: Generator foo has no cost information associated with it']) + + def test_pw_repeat_value_and_cost(self): + curve = example_pw_curve() + orig_values = copy.deepcopy(curve['values']) + curve['values'].insert(2, (30, _example_quadratic(30))) + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, orig_values) + + def test_pw_repeat_value(self): + curve = example_pw_curve() + curve['values'].insert(2, (30, _example_quadratic(40))) + with self.assertRaises(ValueError): + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + + def test_pw_nonconvex(self): + curve = example_pw_curve() + curve['values'].insert(2, (35, _example_quadratic(20))) + with self.assertRaises(ValueError): + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + + def test_pw_low_p_min(self): + curve = example_pw_curve() + expected_values = copy.deepcopy(curve['values']) + expected_values.pop(0) + expected_values.insert(0, (5, 3)) + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=5, + p_max=90, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, expected_values) + self.assertIsNot(cleaned_values, curve['values']) + + def test_pw_high_p_max(self): + curve = example_pw_curve() + expected_values = copy.deepcopy(curve['values']) + expected_values.pop(-1) + expected_values.append((95, 543)) + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=95, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, expected_values) + self.assertIsNot(cleaned_values, curve['values']) + + def test_extra_pw_pieces_below_pmin(self): + curve = example_pw_curve() + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=30, + p_max=90, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, curve['values'][1:]) + self.assertIsNot(cleaned_values, curve['values']) + + def test_extra_pw_pieces_above_pmax(self): + curve = example_pw_curve() + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=70, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, curve['values'][:-1]) + self.assertIsNot(cleaned_values, curve['values']) + + def test_pw_repeated_slope(self): + curve = dict() + curve['data_type'] = 'cost_curve' + curve['cost_curve_type'] = 'piecewise' + curve['values'] = [(10, 20), + (30, 40), + (50, 60), + (70, 80), + (90, 100)] + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=10, + p_max=90, + gen_name='foo', + t=None) + expected_values = [(10, 20), (90, 100)] + self.assertEqual(cleaned_values, expected_values) + self.assertIsNot(expected_values, curve['values']) + + def test_poly_simple(self): + curve = example_poly_curve() + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=15, + p_max=85, + gen_name='foo', + t=None) + self.assertEqual(cleaned_values, curve['values']) + self.assertIsNot(cleaned_values, curve['values']) + + def test_poly_nonconvex(self): + curve = example_poly_curve() + curve['values'][2] = -1 + with self.assertRaises(ValueError): + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=15, + p_max=85, + gen_name='foo', + t=None) + + def test_poly_cubic(self): + curve = example_poly_curve() + curve['values'][3] = 1 + with self.assertRaises(ValueError): + cleaned_values = tx_utils.validate_and_clean_cost_curve(curve=curve, + curve_type='cost_curve', + p_min=15, + p_max=85, + gen_name='foo', + t=None) diff --git a/egret/model_library/transmission/tx_utils.py b/egret/model_library/transmission/tx_utils.py index 63954047..8c464803 100644 --- a/egret/model_library/transmission/tx_utils.py +++ b/egret/model_library/transmission/tx_utils.py @@ -13,7 +13,13 @@ """ import egret.model_library.transmission.tx_calc as tx_calc +import math from math import radians +import logging +import copy + + +logger = logging.getLogger(__name__) def dicts_of_vr_vj(buses): @@ -539,3 +545,100 @@ def _convert_modeldata_pu(model_data, transform_func, inplace): def radians_from_degrees_dict(d): return {k:radians(d[k]) for k in d} + + +def validate_and_clean_cost_curve(curve, curve_type, p_min, p_max, gen_name, t=None): + """ + Parameters + ---------- + curve: dict + curve_type: str + p_min: float + p_max: float + gen_name: str + t: None or int + + Returns + ------- + cleaned_values: list or dict + If the curve is piecewise, this function will return a list. If the curve is + a polynomial, this function will return a dict. + """ + # validate that what we have is a cost_curve + if curve['data_type'] != curve_type: + raise ValueError(f"cost curve must be of data_type {curve_type}.") + + # get the values, check against something empty + values = curve['values'] + if len(values) == 0: + if t is None: + logger.warning(f"WARNING: Generator {gen_name} has no cost information associated with it") + else: + logger.warning(f"WARNING: Generator {gen_name} has no cost information associated with it at time {t}") + return copy.copy(values) + + # if we have a piecewise cost curve, ensure its convexity past p_min + # if no curve_type+'_type' is specified, we assume piecewise (for backwards + # compatibility with no 'fuel_curve_type') + if curve_type + '_type' not in curve or \ + curve[curve_type + '_type'] == 'piecewise': + cleaned_values = list() + last_slope = None + for (o1, c1), (o2, c2) in zip(values, values[1:]): + if o2 <= p_min or math.isclose(p_min, o2): + continue + + if p_max <= o1 or math.isclose(p_max, o1): + continue + + if len(cleaned_values) == 0: + cleaned_values.append((o1, c1)) + + if math.isclose(o2, o1): + if math.isclose(c2, c1): + continue + raise ValueError(f"Found piecewise {curve_type} for generator {gen_name} at time {t} " + + "with nearly infinite slope.") + + cleaned_values.append((o2, c2)) + + # else p_min < o2 + if last_slope is None: + last_slope = (c2 - c1) / (o2 - o1) + continue + this_slope = (c2 - c1) / (o2 - o1) + if this_slope < last_slope and not math.isclose(this_slope, last_slope): + raise ValueError(f"Piecewise {curve_type} must be convex above p_min. " + + f"Found non-convex piecewise {curve_type} for generator {gen_name} at time {t}") + + # combine pieces if slope is the same + if math.isclose(this_slope, last_slope): + cleaned_values.pop(-2) + + # match first and last point with p_min and p_max + first_slope = (cleaned_values[1][1] - cleaned_values[0][1]) / (cleaned_values[1][0] - cleaned_values[0][0]) + first_intercept = cleaned_values[0][1] - first_slope * cleaned_values[0][0] + first_cost = first_slope * p_min + first_intercept + cleaned_values.pop(0) + cleaned_values.insert(0, (p_min, first_cost)) + + last_slope = (cleaned_values[-1][1] - cleaned_values[-2][1]) / (cleaned_values[-1][0] - cleaned_values[-2][0]) + last_intercept = cleaned_values[-1][1] - last_slope * cleaned_values[-1][0] + last_cost = last_slope * p_max + last_intercept + cleaned_values.pop() + cleaned_values.append((p_max, last_cost)) + + # if we have a quadratic cost curve, ensure its convexity + elif curve[curve_type + '_type'] == 'polynomial': + if set(values.keys()) <= {0, 1, 2}: + if 2 in values and values[2] < 0: + raise ValueError(f"Polynomial {curve_type}s must be convex. " + + f"Found non-convex {curve_type} for generator {gen_name} at time {t}.") + else: + raise ValueError(f"Polynomial {curve_type}s must be quadratic. " + + f"Found non-quatric {curve_type} for generator {gen_name} at time {t}.") + cleaned_values = copy.copy(values) + else: + raise Exception(f"Unexpected {curve_type}_type") + + return cleaned_values diff --git a/egret/model_library/unit_commitment/params.py b/egret/model_library/unit_commitment/params.py index 0f1eaa39..9365cd2f 100644 --- a/egret/model_library/unit_commitment/params.py +++ b/egret/model_library/unit_commitment/params.py @@ -889,68 +889,25 @@ def _check_curve(m, g, curve, curve_type): ## first, get a cost_curve out of time series if curve['data_type'] == 'time_series': curve_t = curve['values'][i] + _t = t else: - curve_t = curve - - ## validate that what we have is a cost_curve - if curve_t['data_type'] != curve_type: - raise Exception("p_cost must be of data_type cost_curve.") - - ## get the values, check against something empty - values = curve_t['values'] - if len(values) == 0: - if curve_t == curve: - logger.warning("WARNING: Generator {} has no cost information associated with it".format(g)) - return True - else: - logger.warning("WARNING: Generator {} has no cost information associated with it at time {}".format(g,t)) + curve_t = curve + _t = None - ## if we have a piecewise cost curve, ensure its convexity past p_min - ## if no curve_type+'_type' is specified, we assume piecewise (for backwards - ## compatibility with no 'fuel_curve_type') - if curve_type+'_type' not in curve_t or \ - curve_t[curve_type+'_type'] == 'piecewise': - p_min = value(m.MinimumPowerOutput[g,t]) - last_slope = None - for (o1, c1), (o2, c2) in zip(values, values[1:]): - if o2 <= p_min or math.isclose(p_min, o2): - continue - if math.isclose(o2,o1): - if math.isclose(c2,c1): - continue - raise Exception("Piecewise {} must be convex above p_min. ".format(curve_type) + \ - "Found non-convex piecewise {} for generator {} at time {}".format(curve_type,g,t)) - ## else p_min > o2 - if last_slope is None: - last_slope = (c2-c1)/(o2-o1) - continue - this_slope = (c2-c1)/(o2-o1) - if this_slope < last_slope and not math.isclose(this_slope, last_slope): - raise Exception("Piecewise {} must be convex above p_min. ".format(curve_type) + \ - "Found non-convex piecewise {} for generator {} at time {}".format(curve_type,g,t)) - ## verify the last output value is at least p_max - o_last = values[-1][0] - if value(m.MaximumPowerOutput[g,t]) > o_last and \ - not math.isclose(value(m.MaximumPowerOutput[g,t]), o_last): - raise Exception("{} does not contain p_max for generator {} at time {}".format(curve_type,g,t)) - - ## if we have a quadratic cost curve, ensure its convexity - elif curve_t[curve_type+'_type'] == 'polynomial': + tx_utils.validate_and_clean_cost_curve(curve_t, + curve_type=curve_type, + p_min=value(m.MinimumPowerOutput[g, t]), + p_max=value(m.MaximumPowerOutput[g, t]), + gen_name=g, + t=_t) + + if curve_t[curve_type + '_type'] == 'polynomial': if not _check_curve.warn_piecewise_approx: logger.warning("WARNING: Polynomial cost curves will be approximated using piecewise segments") _check_curve.warn_piecewise_approx = True - values = curve_t['values'] - if set(values.keys()) <= {0,1,2}: - if 2 in values and values[2] < 0: - raise Exception("Polynomial {}s must be convex. ".format(curve_type) + \ - "Found non-convex {} for generator {} at time {}.".format(curve_type,g,t)) - if curve_t == curve: ## in this case, no need to check the other time periods - return - else: - raise Exception("Polynomial {}s must be quatric. ".format(curve_type) + \ - "Found non-quatric {} for generator {} at time {}.".format(curve_type,g,t)) - else: - raise Exception("Unexpected {}_type".format(curve_type)) + + if curve['data_type'] != 'time_series': + break ## set "static" variable for this function _check_curve.warn_piecewise_approx = False diff --git a/egret/models/acopf.py b/egret/models/acopf.py index 530fc69a..ae9cee79 100644 --- a/egret/models/acopf.py +++ b/egret/models/acopf.py @@ -67,7 +67,7 @@ def _validate_and_extract_slack_penalties(model_data): return model_data.data['system']['load_mismatch_cost'], model_data.data['system']['q_load_mismatch_cost'] -def _create_base_power_ac_model(model_data, include_feasibility_slack=False): +def _create_base_power_ac_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'): md = model_data.clone_in_service() tx_utils.scale_ModelData_to_pu(md, inplace=True) @@ -229,28 +229,30 @@ def _create_base_power_ac_model(model_data, include_feasibility_slack=False): ) # declare the generator cost objective - pw_cost_gens = list() - p_cost = gen_attrs['p_cost'] - if 'q_cost' in gen_attrs: - q_cost = gen_attrs['q_cost'] - else: - q_cost = dict() - for gen_name in gen_attrs['names']: - if p_cost[gen_name]['cost_curve_type'] == 'piecewise' or (gen_name in q_cost and q_cost[gen_name]['cost_curve_type'] == 'piecewise'): - pw_cost_gens.append(gen_name) - - if len(pw_cost_gens) > 0: - libgen.declare_var_pg_cost(model=model, index_set=pw_cost_gens) - libgen.declare_var_qg_cost(model=model, index_set=pw_cost_gens) - libgen.declare_piecewise_cost_cons(model=model, index_set=pw_cost_gens, p_costs=p_cost, q_costs=q_cost) - - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'], - q_costs=gen_attrs.get('q_cost', None)) - + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) - obj_expr += sum(model.qg_operating_cost[gen_name] for gen_name in model.qg_operating_cost) + q_costs = gen_attrs.get('q_cost', None) + if q_costs is not None: + pw_qg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=q_costs)) + if len(pw_qg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_qg(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + libgen.declare_qg_delta_qg_con(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + else: + libgen.declare_var_qg_cost(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + libgen.declare_piecewise_qg_cost_cons(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + libgen.declare_expression_qg_operating_cost(model=model, index_set=gen_attrs['names'], q_costs=q_costs, pw_formulation=pw_cost_model) + obj_expr += sum(model.qg_operating_cost[gen_name] for gen_name in model.qg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr diff --git a/egret/models/dcopf.py b/egret/models/dcopf.py index e413f764..adb4c99c 100644 --- a/egret/models/dcopf.py +++ b/egret/models/dcopf.py @@ -49,7 +49,7 @@ def _include_feasibility_slack(model, bus_names, bus_p_loads, gens_by_bus, gen_a for bus_name in bus_names) return p_rhs_kwargs, penalty_expr -def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, include_feasibility_slack=False): +def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, include_feasibility_slack=False, pw_cost_model='delta'): md = model_data.clone_in_service() tx_utils.scale_ModelData_to_pu(md, inplace = True) @@ -191,13 +191,19 @@ def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, inclu coordinate_type=CoordinateType.POLAR ) - ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'] - ) - + # declare the generator cost objective + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr @@ -205,7 +211,7 @@ def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, inclu return model, md -def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_point=BasePointType.FLATSTART, ptdf_options=None): +def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta'): ptdf_options = lpu.populate_default_ptdf_options(ptdf_options) @@ -342,13 +348,19 @@ def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_po approximation_type=ApproximationType.PTDF, ) - ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'] - ) - + # declare the generator cost objective + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr diff --git a/egret/models/tests/test_dcopf.py b/egret/models/tests/test_dcopf.py index 1f678b02..c8ff17f2 100644 --- a/egret/models/tests/test_dcopf.py +++ b/egret/models/tests/test_dcopf.py @@ -17,6 +17,7 @@ from egret.models.dcopf import * from egret.data.model_data import ModelData from parameterized import parameterized +import numpy as np current_dir = os.path.dirname(os.path.abspath(__file__)) case_names = ['pglib_opf_case3_lmbd','pglib_opf_case30_ieee','pglib_opf_case300_ieee'] @@ -72,5 +73,86 @@ def test_ptdf_dcopf_model(self, test_case, soln_case, include_kwargs=False): comparison = math.isclose(md.data['system']['total_cost'], md_soln.data['system']['total_cost'], rel_tol=1e-6) self.assertTrue(comparison) + +def poly_cost_to_pw_cost(md: ModelData, num_points=10): + gen_attrs = md.attributes(element_type='generator') + p_cost = gen_attrs['p_cost'] + for gen_name, cost_dict in p_cost.items(): + assert cost_dict['cost_curve_type'] == 'polynomial' + cost_dict['cost_curve_type'] = 'piecewise' + poly_coefs = cost_dict['values'] + pw_values = list() + p_min = gen_attrs['p_min'][gen_name] + p_max = gen_attrs['p_max'][gen_name] + for pt in np.linspace(p_min, p_max, num_points): + cost = sum(coef*pt**exponent for exponent, coef in poly_coefs.items()) + pw_values.append((pt, cost)) + cost_dict['values'] = pw_values + + +class TestPWCostBTheta(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + download_dir = os.path.join(current_dir, 'transmission_test_instances') + if not os.path.exists(os.path.join(download_dir, 'pglib-opf-master')): + from egret.thirdparty.get_pglib_opf import get_pglib_opf + get_pglib_opf(download_dir) + + def test_case30_pw_cost_b_theta(self): + original_md = ModelData.read(os.path.join(current_dir, 'transmission_test_instances', 'pglib-opf-master', 'pglib_opf_case30_as.m')) + + md = original_md.clone() + m_poly, scaled_md = create_btheta_dcopf_model(md) + opt = pe.SolverFactory('ipopt') + res = opt.solve(m_poly) + poly_obj = pe.value(m_poly.obj.expr) + self.assertAlmostEqual(poly_obj, 767.60209944437236, places=2) + + poly_cost_to_pw_cost(md, num_points=2) + m_pw, scaled_md = create_btheta_dcopf_model(md) + res = opt.solve(m_pw) + pw_obj = pe.value(m_pw.obj.expr) + self.assertAlmostEqual(pw_obj, 788.18275424543515, places=2) + + md = original_md.clone() + poly_cost_to_pw_cost(md, num_points=10) + m_pw, scaled_md = create_btheta_dcopf_model(md) + res = opt.solve(m_pw) + pw_obj = pe.value(m_pw.obj.expr) + self.assertAlmostEqual(pw_obj, 767.74029197558707, places=2) + + +class TestPWCostPTDF(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + download_dir = os.path.join(current_dir, 'transmission_test_instances') + if not os.path.exists(os.path.join(download_dir, 'pglib-opf-master')): + from egret.thirdparty.get_pglib_opf import get_pglib_opf + get_pglib_opf(download_dir) + + def test_case30_pw_cost_b_theta(self): + original_md = ModelData.read(os.path.join(current_dir, 'transmission_test_instances', 'pglib-opf-master', 'pglib_opf_case30_as.m')) + + md = original_md.clone() + m_poly, scaled_md = create_ptdf_dcopf_model(md) + opt = pe.SolverFactory('ipopt') + res = opt.solve(m_poly, tee=False) + poly_obj = pe.value(m_poly.obj.expr) + self.assertAlmostEqual(poly_obj, 767.60209944437236, places=2) + + poly_cost_to_pw_cost(md, num_points=2) + m_pw, scaled_md = create_ptdf_dcopf_model(md) + res = opt.solve(m_pw, tee=False) + pw_obj = pe.value(m_pw.obj.expr) + self.assertAlmostEqual(pw_obj, 783.89649629510222, places=2) + + md = original_md.clone() + poly_cost_to_pw_cost(md, num_points=10) + m_pw, scaled_md = create_ptdf_dcopf_model(md) + res = opt.solve(m_pw, tee=False) + pw_obj = pe.value(m_pw.obj.expr) + self.assertAlmostEqual(pw_obj, 767.74029197558707, places=2) + + if __name__ == '__main__': unittest.main() From 950b550308888303a0d236a509d59e04d8fd34f0 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Wed, 16 Jun 2021 17:06:20 -0600 Subject: [PATCH 4/6] fixing issue with missing 'fuel_curve_type' --- egret/model_library/unit_commitment/params.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/egret/model_library/unit_commitment/params.py b/egret/model_library/unit_commitment/params.py index 9365cd2f..6925b484 100644 --- a/egret/model_library/unit_commitment/params.py +++ b/egret/model_library/unit_commitment/params.py @@ -901,7 +901,10 @@ def _check_curve(m, g, curve, curve_type): gen_name=g, t=_t) - if curve_t[curve_type + '_type'] == 'polynomial': + # if no curve_type+'_type' is specified, we assume piecewise (for backwards + # compatibility with no 'fuel_curve_type') + if curve_type + '_type' in curve and \ + curve_t[curve_type + '_type'] == 'polynomial': if not _check_curve.warn_piecewise_approx: logger.warning("WARNING: Polynomial cost curves will be approximated using piecewise segments") _check_curve.warn_piecewise_approx = True From 714ac0544f488b4500a4c36591a9c674c7bddf75 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Wed, 16 Jun 2021 18:40:42 -0600 Subject: [PATCH 5/6] removing declare_expression_pgqg_operating_cost --- egret/models/acopf.py | 31 ++++++++++++++++++++------- egret/models/copperplate_dispatch.py | 16 +++++++++----- egret/models/dcopf_losses.py | 32 +++++++++++++++++++--------- egret/models/scopf.py | 16 +++++++++----- 4 files changed, 67 insertions(+), 28 deletions(-) diff --git a/egret/models/acopf.py b/egret/models/acopf.py index ae9cee79..85bdb00d 100644 --- a/egret/models/acopf.py +++ b/egret/models/acopf.py @@ -599,17 +599,32 @@ def create_riv_acopf_model(model_data, include_feasibility_slack=False): coordinate_type=CoordinateType.RECTANGULAR) ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'], - q_costs=gen_attrs.get('q_cost', None) - ) - + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + q_costs = gen_attrs.get('q_cost', None) + if q_costs is not None: + pw_qg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=q_costs)) + if len(pw_qg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_qg(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + libgen.declare_qg_delta_qg_con(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + else: + libgen.declare_var_qg_cost(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + libgen.declare_piecewise_qg_cost_cons(model=model, index_set=pw_qg_cost_gens, q_costs=q_costs) + libgen.declare_expression_qg_operating_cost(model=model, index_set=gen_attrs['names'], q_costs=q_costs, pw_formulation=pw_cost_model) + obj_expr += sum(model.qg_operating_cost[gen_name] for gen_name in model.qg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr - if hasattr(model, 'qg_operating_cost'): - obj_expr += sum(model.qg_operating_cost[gen_name] for gen_name in model.qg_operating_cost) model.obj = pe.Objective(expr=obj_expr) diff --git a/egret/models/copperplate_dispatch.py b/egret/models/copperplate_dispatch.py index d5d480d8..38645e02 100644 --- a/egret/models/copperplate_dispatch.py +++ b/egret/models/copperplate_dispatch.py @@ -94,12 +94,18 @@ def create_copperplate_dispatch_approx_model(model_data, include_feasibility_sla ) ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'] - ) - + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr diff --git a/egret/models/dcopf_losses.py b/egret/models/dcopf_losses.py index 15a6c0fb..c6876ca4 100644 --- a/egret/models/dcopf_losses.py +++ b/egret/models/dcopf_losses.py @@ -174,12 +174,18 @@ def create_btheta_losses_dcopf_model(model_data, relaxation_type=RelaxationType. ) ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'] - ) - + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr @@ -304,12 +310,18 @@ def create_ptdf_losses_dcopf_model(model_data, include_feasibility_slack=False, ) ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'] - ) - + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr diff --git a/egret/models/scopf.py b/egret/models/scopf.py index 15dac032..883c4f73 100644 --- a/egret/models/scopf.py +++ b/egret/models/scopf.py @@ -168,12 +168,18 @@ def create_scopf_model(model_data, include_feasibility_slack=False, base_point=B lpu.add_initial_monitored_branches(model, branches, branches_idx, ptdf_options, PTDF) ### declare the generator cost objective - libgen.declare_expression_pgqg_operating_cost(model=model, - index_set=gen_attrs['names'], - p_costs=gen_attrs['p_cost'] - ) - + p_costs = gen_attrs['p_cost'] + pw_pg_cost_gens = list(libgen.pw_gen_generator(gen_attrs['names'], costs=p_costs)) + if len(pw_pg_cost_gens) > 0: + if pw_cost_model == 'delta': + libgen.declare_var_delta_pg(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_pg_delta_pg_con(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + else: + libgen.declare_var_pg_cost(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_piecewise_pg_cost_cons(model=model, index_set=pw_pg_cost_gens, p_costs=p_costs) + libgen.declare_expression_pg_operating_cost(model=model, index_set=gen_attrs['names'], p_costs=p_costs, pw_formulation=pw_cost_model) obj_expr = sum(model.pg_operating_cost[gen_name] for gen_name in model.pg_operating_cost) + if include_feasibility_slack: obj_expr += penalty_expr From 4997d573f4f8333e15bc123c1a7092172e36c5e7 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Thu, 17 Jun 2021 08:08:24 -0600 Subject: [PATCH 6/6] making pw_cost_model a top-level argument --- egret/models/acopf.py | 17 ++++++++++------- egret/models/copperplate_dispatch.py | 2 +- egret/models/dcopf_losses.py | 6 ++++-- egret/models/scopf.py | 3 ++- 4 files changed, 17 insertions(+), 11 deletions(-) diff --git a/egret/models/acopf.py b/egret/models/acopf.py index 85bdb00d..7e793bdb 100644 --- a/egret/models/acopf.py +++ b/egret/models/acopf.py @@ -261,8 +261,9 @@ def _create_base_power_ac_model(model_data, include_feasibility_slack=False, pw_ return model, md -def create_atan_acopf_model(model_data, include_feasibility_slack=False): - model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack) +def create_atan_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'): + model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack, + pw_cost_model=pw_cost_model) branch_attrs = md.attributes(element_type='branch') bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus']) @@ -304,8 +305,9 @@ def create_atan_acopf_model(model_data, include_feasibility_slack=False): return model, md -def create_psv_acopf_model(model_data, include_feasibility_slack=False): - model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack) +def create_psv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'): + model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack, + pw_cost_model=pw_cost_model) bus_attrs = md.attributes(element_type='bus') branch_attrs = md.attributes(element_type='branch') bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus']) @@ -348,8 +350,9 @@ def create_psv_acopf_model(model_data, include_feasibility_slack=False): return model, md -def create_rsv_acopf_model(model_data, include_feasibility_slack=False): - model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack) +def create_rsv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'): + model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack, + pw_cost_model=pw_cost_model) bus_attrs = md.attributes(element_type='bus') branch_attrs = md.attributes(element_type='branch') bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus']) @@ -390,7 +393,7 @@ def create_rsv_acopf_model(model_data, include_feasibility_slack=False): return model, md -def create_riv_acopf_model(model_data, include_feasibility_slack=False): +def create_riv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'): md = model_data.clone_in_service() tx_utils.scale_ModelData_to_pu(md, inplace = True) diff --git a/egret/models/copperplate_dispatch.py b/egret/models/copperplate_dispatch.py index 38645e02..1c36dc08 100644 --- a/egret/models/copperplate_dispatch.py +++ b/egret/models/copperplate_dispatch.py @@ -44,7 +44,7 @@ def _validate_and_extract_slack_penalty(model_data): assert('load_mismatch_cost' in model_data.data['system']) return model_data.data['system']['load_mismatch_cost'] -def create_copperplate_dispatch_approx_model(model_data, include_feasibility_slack=False): +def create_copperplate_dispatch_approx_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'): md = model_data.clone_in_service() tx_utils.scale_ModelData_to_pu(md, inplace = True) diff --git a/egret/models/dcopf_losses.py b/egret/models/dcopf_losses.py index c6876ca4..6ea7c224 100644 --- a/egret/models/dcopf_losses.py +++ b/egret/models/dcopf_losses.py @@ -33,7 +33,8 @@ from math import pi, radians, degrees -def create_btheta_losses_dcopf_model(model_data, relaxation_type=RelaxationType.SOC, include_angle_diff_limits=False, include_feasibility_slack=False): +def create_btheta_losses_dcopf_model(model_data, relaxation_type=RelaxationType.SOC, + include_angle_diff_limits=False, include_feasibility_slack=False, pw_cost_model='delta'): md = model_data.clone_in_service() tx_utils.scale_ModelData_to_pu(md, inplace = True) @@ -194,7 +195,8 @@ def create_btheta_losses_dcopf_model(model_data, relaxation_type=RelaxationType. return model, md -def create_ptdf_losses_dcopf_model(model_data, include_feasibility_slack=False, ptdf_options=None): +def create_ptdf_losses_dcopf_model(model_data, include_feasibility_slack=False, + ptdf_options=None, pw_cost_model='delta'): ptdf_options = lpu.populate_default_ptdf_options(ptdf_options) diff --git a/egret/models/scopf.py b/egret/models/scopf.py index 883c4f73..18efeb76 100644 --- a/egret/models/scopf.py +++ b/egret/models/scopf.py @@ -31,7 +31,8 @@ from math import pi, radians, degrees -def create_scopf_model(model_data, include_feasibility_slack=False, base_point=BasePointType.FLATSTART, ptdf_options=None): +def create_scopf_model(model_data, include_feasibility_slack=False, + base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta'): ptdf_options = lpu.populate_default_ptdf_options(ptdf_options)