4. Optimization Problems¶
This section describes how to formulate power network optimization problems using PFNET.
4.1. Objective Function¶
The objective function \(\phi\) for a network optimization problem created using PFNET is of the form
where \(w_i\) are weights, \(\varphi_i\) are general linear or nonlinear functions, and \(x\) is a vector of network variables. Each weight-function pair in the summation is represented by an object of type FunctionBase. An instance of this type can be constructed from the Function class, which requires specifying the function name, weight, and the Network object to be associated with the function. The following example sets all bus voltage magnitudes as variables and constructs a function that penalizes voltage magnitude deviations from ideal values:
>>> import pfnet
>>> net = net = pfnet.PyParserMAT().parse('ieee14.m')
>>> net.set_flags('bus',
... 'variable',
... 'any',
... 'voltage magnitude')
>>> func = pfnet.Function('voltage magnitude regularization', 0.3, net)
>>> print func.name == 'voltage magnitude regularization'
True
>>> print func.weight
0.3
After a FunctionBase object is created, its value, gradient, and Hessian are zero, an empty vector, and an empty matrix, respectively. Before evaluating the function at a specific vector of values, it must be analyzed using the FunctionBase class method analyze(). This routine analyzes the function and allocates the required vectors and matrices for storing its gradient and Hessian. After this, the function can be evaluated using the method eval():
>>> x = net.get_var_values()
>>> func.analyze()
>>> func.eval(x)
The value \(\varphi_i(x)\), gradient \(\nabla \varphi_i(x)\), and Hessian \(\nabla^2 \varphi_i(x)\) of a function can then be extracted from the phi, gphi and Hphi attributes, respectively:
>>> print x.shape
(14,)
>>> print func.phi
0.255
>>> print type(func.gphi), func.gphi.shape
<type 'numpy.ndarray'> (14,)
>>> print type(func.Hphi), func.Hphi.shape
<class 'scipy.sparse.coo.coo_matrix'> (14, 14)
For the Hessian matrix, only the lower triangular part is stored.
Details about each of the different functions available in PFNET are provided below.
4.1.1. Active power generation cost¶
This function is associated with the string "generation cost". It measures active power generation cost by the expression
where \(P_k(t)\) are generator active powers in per unit system base power, \(t\) are time periods, and \(q_{k0}\), \(q_{k1}\), and \(q_{k2}\) are constant coefficients. These coefficients correspond to the attributes cost_coeff_Q0, cost_coeff_Q1, and cost_coeff_Q2 of each Generator object.
4.1.2. Active power consumption utility¶
This function is associated with the string "consumption utility". It measures active power consumption utility by the expression
where \(P_k(t)\) are load active powers in per unit system base power, \(t\) are time periods, and \(q_{k0}\), \(q_{k1}\), and \(q_{k2}\) are constant coefficients. These coefficients correspond to the attributes util_coeff_Q0, util_coeff_Q1, and util_coeff_Q2 of each Load object.
4.1.3. Net Active Power Consumption Cost¶
This function is associated with the string "net consumption cost". It measures the total cost of net active power consumption over the time periods using the price defined by the price attribute of each Bus object.
4.1.4. Voltage magnitude regularization¶
This function is associated with the string "voltage magnitude regularization". It penalizes deviations of bus voltage magnitudes from ideal values. It is defined by the expression
where \(t\) are time periods, \(v\) are bus voltage magnitudes, \(v^s\) are voltage magnitude set points (one for buses not regulated by generators), and \(\Delta v\) is a normalization factor.
4.1.5. Voltage angle regularization¶
This function is associated with the string "voltage angle regularization". It penalizes large bus voltage angles and voltage angle differences across branches. It is defined by the expression
where \(t\) are time periods, \(\theta\) are bus voltage angles, \(\phi\) are branch phase shifts, and \(\Delta \theta\) is a normalization factor. Only terms that include optimization variables are included in the summation.
4.1.6. Generator powers regularization¶
This function is associated with the string "generator powers regularization". It penalizes deviations of generator powers from the midpoint of their ranges. It is defined by the expression
where \(t\) are time periods, \(P^g\) and \(Q^g\) are generator active and reactive powers, \(\bar{P}\) and \(\bar{Q}\) are midpoints of generator active and reactive power ranges, and \(\Delta P = \Delta Q\) are normalization factors. Only terms that include optimization variables are included in the summation.
4.1.7. Transformer tap ratio regularization¶
This function is associated with the string "tap ratio regularization". It penalizes deviations of tap ratios of tap-changing transformers from their initial values. It is defined by the expression
where \(\tau\) are time periods, \(t\) are tap ratios of tap-changing transformers, \(t^0\) are their initial values, and \(\Delta t\) is a normalization factor.
4.1.8. Transformer phase shift regularization¶
This function is associated with the string "phase shift regularization". It penalizes deviations of phase shifts of phase shifting transformers from their initial value. It is defined by the expression
where \(t\) are time periods, \(\phi\) are phase shifts of phase-shifting transformers, \(\phi^0\) are their initial values, and \(\Delta \phi\) is a normalization factor. Only terms that include optimization variables are included in the summation.
4.1.9. Switched shunt susceptance regularization¶
This function is associated with the string "susceptance regularization". It penalizes deviations of susceptances of switched shunt devices from their initial values. It is defined by the expression
where \(t\) are time periods, \(b\) are susceptances of switched shunt devices, \(b^0\) are their initial values, and \(\Delta b\) is a normalization factor. Only terms that include optimization variables are included in the summation.
4.1.10. Voltage magnitude soft limit penalty¶
This function is associated with the string "soft voltage magnitude limits". It penalizes deviations of bus voltage magnitudes from the mid point of their ranges. It is defined by the expression
where \(t\) are time periods, \(v\) are bus voltage magnitudes, \(\bar{v}\) are the mid points of their ranges, and \(\Delta v\) is a normalization factor. Only terms that include optimization variables are included in the summation.
4.1.11. Sparsity inducing penalty for controls¶
This function is associated with the string "sparse controls penalty". It encourages sparse control adjustments with the expression
where \(t\) are time periods, \(u\) are control quantities, \(u^0\) are their current values, and \(\epsilon\) is a small positive scalar. The normalization factors \(\Delta u_k\) are given by
where \(u^{\max}\) and \(u^{\min}\) are control limits, and \(\delta\) is a small positive scalar. The control quantities that are considered by this function are specified using the Network class methods set_flags() or set_flags_of_component() using the flag "sparse".
4.1.12. General Variable Regularization¶
This function is associated with the string "variable regularization". It consists of the quadratic function
where the vectors \(x_0\) and \(w\) are parameters that can be set using the FunctionBase class method set_parameter(), as shown below:
>>> import numpy as np
>>> func = pfnet.Function('variable regularization', 1., net)
>>> w = np.random.randn(net.num_vars)
>>> x0 = np.random.randn(net.num_vars)
>>> func.set_parameter('w', w)
>>> func.set_parameter('x0', x0)
4.2. Constraints¶
Constraints in PFNET are of the form
where \(A\) and \(G\) are sparse matrices, \(b\), \(l\) and \(u\) are vectors, \(f\) is a vector-valued nonlinear function, and \(x\) and \(y\) are vectors of network variables and extra or auxiliary constraint variables, respectively. They are represented by objects of type ConstraintBase. An instance of this type can be constructed from the class Constraint, whose constructor requires specifying the constraint name and the Network to be associated with the constraint. The following example sets all bus voltage magnitudes and angles as variables and constructs the AC power balance constraints:
>>> import pfnet
>>> pfnet.PyParserMAT().parse('ieee14.m')
>>> net.set_flags('bus',
... 'variable',
... 'any',
... ['voltage magnitude','voltage angle'])
>>> print net.num_vars == 2*net.num_buses
True
>>> constr = pfnet.Constraint('AC power balance',net)
>>> print constr.name == 'AC power balance'
True
Before a ConstraintBase object can be used, it must be initialized using the ConstraintBase class method analyze(). This routine analyzes the constraint and allocates the required vectors and matrices. After this, the number of extra variables can be obtained from the attribute num_extra_vars, and the constraint can be evaluated using the method eval():
>>> x = net.get_var_values()
>>> constr.analyze()
>>> print constr.num_extra_vars
0
>>> constr.eval(x)
The matrices and vectors associated with the linear constraints can be extracted from the A, G, b, l and u attributes of the ConstraintBase object. The vector of violations and Jacobian matrix of the nonlinear constraints can be extracted from the attributes f and J, respectively. Also, the Hessian matrix of any individual nonlinear constraint function \(f_i(x)\) can be extracted using the class method get_H_single(). The following example shows how to extract the largest power flow mismatch in per unit system base power and the Hessian matrix corresponding to the active power balance constraint of a bus:
>>> import numpy as np
>>> f = constr.f
>>> print type(f), f.shape
<type 'numpy.ndarray'> (28,)
>>> print np.linalg.norm(f,np.inf)
0.042
>>> bus = net.get_bus(5)
>>> Hi = constr.get_H_single(bus.dP_index)
>>> print type(Hi), Hi.shape, Hi.nnz
<class 'scipy.sparse.coo.coo_matrix'> (28, 28) 27
As before, all Hessian matrices have stored only the lower triangular part. In addition to being possible to extract Hessian matrices of individual nonlinear constraints, it is also possible to construct any linear combination of these individual Hessian matrices. This can be done using the ConstraintBase class method combine_H(). After this, the resulting matrix can be extracted from the H_combined attribute:
>>> coefficients = np.random.randn(f.size)
>>> constr.combine_H(coefficients)
>>> H = constr.H_combined
>>> print type(H), H.shape, H.nnz
<class 'scipy.sparse.coo.coo_matrix'> (28, 28) 564
Lagrange multiplier estimates of the linear and nonlinear constraints can be used to store sensitivity information in the network components associated with the constraints. This is done using the class method store_sensitivities(). Component-specific attributes that store sensitivity information are described in the API Reference section.
Lastly, the ConstraintBase class provides methods for extracting information about the rows of each linear and nonlinear constraint. These methods are get_A_row_info_string(), get_J_row_info_string(), and get_G_row_info_string().
Details about each of the different constraints available in PFNET are provided below.
4.2.1. AC Power balance¶
This constraint is associated with the string "AC power balance". It enforces “AC” active and reactive power balance at every bus of the network. It is given by
where \(t\) are time periods, \(P^g\) and \(Q^g\) are generator active and reactive powers, \(P^l\) and \(Q^l\) are load active and reactive powers, \(S^{sh}\) are apparent powers flowing out of buses through shunt devices, \(S\) are apparent powers flowing out of buses through branches, \(n\) is the number of buses, \(T\) is the number of time periods, and \([n] := \{1,\ldots,n\}\).
The Bus class attributes dP_index and dQ_index can be used to obtain the row indices of the constraint function f and Jacobian J that are associated with the active and reactive power mismatches of a specific bus.
4.2.2. DC Power balance¶
This constraint is associated with the string "DC power balance". It enforces “DC” active power balance at every bus of the network. It is given by
where \(t\) are time periods, \(P^g\) are generator active powers, \(P^l\) are load active powers, \(b_{km}\) are branch susceptances, \(\theta_k\) are bus voltage angles, \(\phi_{km}\) are phase shifts of phase-shifting transformers, \(n\) is the number of buses, \(T\) is the number of time periods, and \([n] := \{1,\ldots,n\}\).
The Bus class attribute dP_index can be used to obtain the row indices of the constraint matrix A and right-hand-side b that are associated with the active (DC) power mismatches of a specific bus.
4.2.3. Linearized AC Power balance¶
This constraint is associated with the string "linearized AC power balance". It enforces active and reactive power balance at every bus of the network using a first-order Taylor expansion of the AC power balance constraints. It is given by
where \(x_0\) is the vector of current variable values, \(f(x_0)\) is the vector of AC bus power mismatches, and \(J(x_0)\) is the Jacobian of \(f\) at \(x_0\).
4.2.4. DC branch flow limits¶
This constraint is associated with the string "DC branch flow limits". It enforces branch “DC” power flow limits due to thermal ratings. It is given by
for each branch \((k,m)\) and time period \(t\), where \(b_{km}\) are branch susceptances, \(\theta_k\) are bus voltage angles, \(\phi_{km}\) are phase shifts of phase-shifting transformers, and \(P^{\max}_{km}\) are branch power flow limits.
4.2.5. AC branch flow limits¶
This constraint is associated with the string "AC branch flow limits". It enforces branch “AC” power flow limits due to thermal ratings based on current magnitudes. It is given by
for each time period \(t\) and branch \((k,m)\) with \(i^{\max}_{km} \neq 0\) (thermal rating A), where \(i_{km}\) denotes the current leaving bus \(k\) towards bus \(m\). This constraint utilizes auxiliary variables (slacks) in order to be expressed in the standard form described at the beginning of Section Constraints.
4.2.6. Linearized AC branch flow limits¶
This constraint is associated with the string "linearized AC branch flow limits". It enforces branch “AC” power flow limits due to thermal ratings based on current magnitudes using conservative linear constraints constructed with the external library Line Flow [DS2017] using default parameters.
4.2.7. Variable fixing¶
This constraint is associated with the string "variable fixing". It constrains specific variables to be fixed at their current value. The variables to be fixed are specified using the Network class methods set_flags() or set_flags_of_component() with the flag "fixed".
4.2.8. Variable bounds¶
This constraint is associated with the string "variable bounds". It constrains specific variables to be inside their bounds. The variables to be bounded are specified using the Network class methods set_flags() or set_flags_of_component() with the flag "bounded".
4.2.9. Generator participation¶
These constraint are associated with the strings "generator active power participation" and "PVPQ switching". They enforce specific active power participations among slack generators connected to the same bus, and reactive power participations among generators regulating the same bus voltage magnitude, respectively. For slack generators, all participate with equal active powers. For voltage-regulating generators, each one participates with a fraction of the combined resources used to regulate the voltage equal to the Generator class field Q_par. The "PVPQ switching" constraint also fixes regulated bus voltage magnitudes to their set points and can modify these constraints depending on whether generator reactive powers reach their limits.
4.2.10. Voltage set-point regulation by generators¶
This constraint is associated with the string "voltage regulation by generators". It enforces voltage set-point regulation by generators. It approximates the constraints
for each time period \(t\) and bus \(k\) whose voltage is regulated by generators, where \(v\) are bus voltage magnitudes, \(v^s\) are their set points, \(v^y\) and \(v^z\) are positive and negative deviations of \(v\) from \(v^s\) (auxiliary variables of the constraint), and \(Q\), \(Q^{\max}\) and \(Q^{\min}\) are aggregate reactive powers and limits of the generators regulating the same bus voltage magnitude.
4.2.11. Voltage band regulation by transformers¶
This constraint is associated with the string "voltage regulation by transformers". It enforces voltage band regulation by tap-changing transformers. It approximates the constraints
for each time period \(\tau\) and bus \(k\) whose voltage is regulated by tap-changing transformers, where \(v\) are bus voltage magnitudes, \(v^{\max}\) and \(v^{\min}\) are their band limits, \(v^l\) and \(v^h\) are voltage violations of band lower and upper limits (auxiliary variables), \(t\) are transformer tap ratios, \(t^0\), \(t^{\max}\) and \(t^{\min}\) are their current values and limits, and \(t^y\) and \(t^z\) are positive and negative deviations of \(t\) from \(t^0\) (auxiliary variables). The above equations assume that the sensitivity between voltage magnitude and transformer tap ratio is positive. If it is negative, \(t^y\) and \(t^z\) are interchanged in the first two complementarity constraints, and \(v^l\) and \(v^h\) are interchanged in the bottom two complementarity constraints.
4.2.12. Voltage band regulation by switched shunts¶
This constraint is associated with the string "voltage regulation by shunts". It enforces voltage band regulation by switched shunt devices. It approximates the constraints
for each time period \(t\) and bus \(k\) whose voltage is regulated by switched shunt devices, where \(v\) are bus voltage magnitudes, \(v^{\max}\) and \(v^{\min}\) are their band limits, \(v^l\) and \(v^h\) are voltage violations of band lower and upper limits (auxiliary variables), \(b\) are switched shunt susceptances, \(b^0\), \(b^{\max}\) and \(b^{\min}\) are their current values and limits, and \(b^y\) and \(b^z\) are positive and negative deviations of \(b\) from \(b^0\) (auxiliary variables).
4.2.13. Generator active power ramp limits¶
This constraint is associated with the string "generator ramp limits". It enforces generator active power ramping limits. It is given by
for each generator \(k\) and time period \(t \in \{1,\ldots,T\}\), where \(P^g\) are generator active powers, and \(\delta P^{\max}\) are generator ramping limits. The ramping limits are defined by the dP_max attribute of each Generator object. For \(t = 1\), \(P^g_k(t-1)\) is the P_prev attribute of a Generator.
4.2.14. Battery dynamics¶
This constraint is associated with the string "battery dynamics". It enforces the dynamic equations of the batteries’ energy levels. It is given by
for each battery \(k\), where \(E\), \(P^c\), and \(P^d\) denote battery energy levels and charging and discharging powers, respectively, and \(E^i\), \(E^f\), \(\eta^c\), and \(\eta^d\) correspond to the attributes E_init , E_final, eta_c, and eta_d of a Battery, respectively. It is noted here that the units of the charging/discharging powers are p.u. system base power, and the units of the energy levels are p.u. system base power times the duration of a time period.
4.2.15. Load Constant Power Factor¶
This constraint is associated with the string "load constant power factor". It forces loads to have a constant power factor at each time. It is given by
for each load \(k\) and time period \(t \in \{1,\ldots,T\}\), where \(P^l\) are active powers, \(Q^l\) are reactive powers, and \(\gamma\) is the load’s target_power_factor.
4.3. Heuristics¶
In power networks, heuristics are often used to model certain controls such as voltage regulation. In PFNET, heuristics are objects of type HeuristicBase. They can be constructed from the Heuristic class, which requires specifying a valid name and a Network object. Heuristics act on a list of constraints and modify the matrices and vectors in order to achieve a desired result. For example, the "PVPQ switching" heuristic modifies the "PVPQ switching" constraint in order to fix a generator reactive power to its limit instead of the regulated voltage to its set point.
4.4. Problems¶
Optimization problems constructed with PFNET are of the form
As already noted, the objective function \(\varphi\) is a weighted sum of functions \(\varphi_i\). The linear and nonlinear constraints \(Ax = b\), \(l \le Gx \le u\), and \(f(x) = 0\) correspond to one or more of the constraints described above. An optimization problem in PFNET is represented by an object of type Problem.
After instantiation, a Problem is empty and one needs to specify the Constraints to include, and the Functions that form the objective function. This can be done using the Problem class methods add_constraint(), and add_function(), respectively. Heuristics can be added using the method add_heuristic(). The following example shows how to construct a simple power flow problem and solve it using the Newton-Raphson method:
import pfnet
from numpy import hstack
from numpy.linalg import norm
from scipy.sparse import bmat
from scipy.sparse.linalg import spsolve
def NRsolve(net):
net.clear_flags()
# bus voltage angles
net.set_flags('bus',
'variable',
'not slack',
'voltage angle')
# bus voltage magnitudes
net.set_flags('bus',
'variable',
'not regulated by generator',
'voltage magnitude')
# slack gens active powers
net.set_flags('generator',
'variable',
'slack',
'active power')
# regulator gens reactive powers
net.set_flags('generator',
'variable',
'regulator',
'reactive power')
p = pfnet.Problem(net)
p.add_constraint(pfnet.Constraint('AC power balance', net))
p.add_constraint(pfnet.Constraint('generator active power participation', net))
p.add_constraint(pfnet.Constraint('PVPQ switching', net))
p.add_heuristic(pfnet.Heuristic('PVPQ switching', net))
p.analyze()
x = p.get_init_point()
p.eval(x)
residual = lambda x: hstack((p.A*x-p.b, p.f))
while norm(residual(x)) > 1e-4:
p.apply_heuristics(x)
x = x + spsolve(bmat([[p.A],[p.J]],format='csr'), -residual(x))
p.eval(x)
net.set_var_values(x)
net.update_properties()
The above routine can then be used as follows:
>>> pfnet.ParserMAT().parse('case3012wp.mat')
>>> print net.bus_P_mis, net.bus_Q_mis
2.79e+0 1.56e+1
>>> NRsolve(net)
>>> print net.bus_P_mis, net.bus_Q_mis
2.37e-6 3.58e-6
As shown in the example, the Problem class method analyze() needs to be called before the vectors and matrices associated with the problem constraints and functions can be used. The method eval() can then be used for evaluating the problem objective and constraint functions at different points. As is the case for Constraints, a Problem has a method combine_H() for forming linear combinations of individual constraint Hessians, and a method store_sensitivities() for storing sensitivity information in the network components associated with the constraints.