diff --git a/egret/model_library/transmission/gen.py b/egret/model_library/transmission/gen.py index fbd0ce05..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,33 +31,236 @@ def declare_var_qg(model, index_set, **kwargs): decl.declare_var('qg', model=model, index_set=index_set, **kwargs) -def declare_expression_pgqg_operating_cost(model, index_set, - p_costs, q_costs=None): +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, 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': + 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) + intercept = cost2 - slope * pt2 + pw_cost_set.add((gen_name, ndx)) + indexed_pw_cost_con[gen_name, ndx] = cost_var >= slope * gen_var + intercept + else: + raise ValueError(f"Unrecognized cost_cureve_type: {cost_dict['cost_curve_type']}") + + +def declare_piecewise_pg_cost_cons(model, index_set, p_costs): + m = model + + m.pg_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 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) + + for gen_name in expr_set: + 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': + 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 + + +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) - 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 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 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': + 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: + m.qg_operating_cost[gen_name] = 0 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..6925b484 100644 --- a/egret/model_library/unit_commitment/params.py +++ b/egret/model_library/unit_commitment/params.py @@ -889,68 +889,28 @@ 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)) - - ## 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': + curve_t = curve + _t = None + + 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 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 - 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 ca16f464..7e793bdb 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) @@ -228,26 +228,42 @@ def _create_base_power_ac_model(model_data, include_feasibility_slack=False): branches=branches ) - ### 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) - ) - + # 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) + 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) 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']) @@ -289,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']) @@ -333,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']) @@ -375,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) @@ -584,17 +602,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..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) @@ -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.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/dcopf_losses.py b/egret/models/dcopf_losses.py index 15a6c0fb..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) @@ -174,12 +175,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 @@ -188,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) @@ -304,12 +312,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..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) @@ -168,12 +169,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 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() 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()