Dynamic Optimal Power Flow in Microgrids using the Alternating Direction Method of Multipliers
Abstract
Smart devices, storage and other distributed technologies have the potential to greatly improve the utilisation of network infrastructure and renewable generation. Decentralised control of these technologies overcomes many scalability and privacy concerns, but in general still requires the underlying problem to be convex in order to guarantee convergence to a global optimum. Considering that AC power flows are nonconvex in nature, and the operation of household devices often requires discrete decisions, there has been uncertainty surrounding the use of distributed methods in a realistic setting. This paper extends prior work on the alternating direction method of multipliers (ADMM) for solving the dynamic optimal power flow (DOPF) problem. We utilise more realistic line and load models, and introduce a twostage approach to managing discrete decisions and uncertainty. Our experiments on a suburbsized microgrid show that this approach provides near optimal results, in a time that is fast enough for receding horizon control. This work brings distributed control of smartgrid technologies closer to reality.
I Introduction
The role of electricity market operators is to supply low cost and reliable electricity to customers. This typically involves solving unit commitment (UC) and/or optimal power flow (OPF) problems at regular intervals. In the classical form of these problems, loads are assumed to be inflexible, so only large generators and network components need to be considered. However, this assumption will no longer hold true when smart devices, distributed generation, storage, and electric vehicles (EVs) are adopted on mass. The traditional centralised market approaches were not designed to either operate where every consumer is an active participant, or handle their unique timecoupled behaviours.
In order to manage this huge increase in problem size, several works [1, 2, 3] have adopted distributed solving techniques. In addition to greatly parallelising the problem, these distributed algorithms help to preserve the privacy of participants. These algorithms require the problem to be convex in order to guarantee convergence to a globally optimal solution. However, the behaviour of many loads at a household level are discrete in nature [4], and the equations that govern how power physically flows on the network are nonconvex.
This paper considers the task of balancing power within a microgrid. We formulate the problem as a dynamic optimal power flow (DOPF) problem to account for multiple time steps over a day, and solve it in a distributed manner by adopting the alternating direction method of multipliers (ADMM) approach presented in [1]. We extend the method to more accurate power flow models and introduce a twostage pricing mechanism to manage integer variables and uncertainty. We find that this approach achieves near optimal solutions in a time frame that is fast enough to work with receding horizon control. This brings the overall approach closer to the point where it can be deployed on a real network.
Ii Related Work
Demand response (DR) is the name often given to the control of distributed technologies. Much of the existing work on DR has focused on using realtime pricing (RTP) as a control signal [5, 6, 7, 8]. In these methods, participants receive a RTP signal and individually optimise their own behaviour, so as to minimise a combination of monetary and discomfort costs. Other approaches have utilised nonpricing control signals, which are simpler to implement, but are limited in the types of loads that they can model [9, 10].
These approaches can be thought of as a form of open loop control, because the agent that sets the control signal (RTP or otherwise), at best can only estimate how consumers will respond to it. The work in [2] takes the idea of RTP and closes the loop. The RTP is iteratively updated by a central agent, with consumers communicating their best responses to the price before implementing them. An alternative iterative procedure is introduced in [3], where consumers cooperate to reduce total generation costs in a distributed manner. This approach is robust to gaming, but individual households only have a small incentive to change their habits, as their costs are largely independent of their own behaviour.
The approaches discussed so far do not model the electricity network, so cannot account for real power losses, reactive power, voltage limits and line thermal limits. Most of the work on distributed algorithms that do model the network have used ADMM as a solving technique, due to its ease in decomposition, and its convergence guarantees on a wide range of problems [11]. One of the first authors to apply ADMM to power networks was Kim et al. [12]. They decomposed an OPF problem into regions and compared ADMM to two other approaches. They found that ADMM provided significant speed improvements over a centralised approach, and that privacy is preserved between regions. In [13], Phan et al. used ADMM to improve the parallel solve time for a securityconstrained OPF (SCOPF) problem. They decomposed across scenarios and found it to have comparable solve times to a Bender’s cut approach in many instances.
The decomposition power of ADMM was taken one step further when Kraning et al. [1] decomposed all network components for a DOPF problem. This is effectively a principled method for settling RTPs for each bus, also known as locational marginal prices. Their experiments showed that very large problems could be solved efficiently in a parallel environment. In this paper we extend this work to include more sophisticated line models, and a method for managing discrete variables and uncertainty, thereby making ADMM a practical approach to modern DOPF problems.
Iii Network Model
The overall objective is to minimise the cost of supplying electricity. We formulate this as DOPF problems embedded within a receding horizon control process, in order to accurately control timecoupled components in an uncertain environment. In receding horizon control, a DOPF is first solved over a horizon of time steps, the decision in the first step is acted on, and then the process repeats with the window shifted forward by one. In this paper we focus on solving the DOPF within a single horizon, and the actions that agents take to implement the first decision. We formalise the network model with this problem in mind.
Note that the notation and formulation that we use is different from what is standard in power systems. This is necessary in order to decompose the network and distribute the problem. Fig. 1 highlights the difference between a standard line diagram and our formulation.
A network consists of a set of components , terminals and connections . Each component (e.g., bus, line, generator, load) has a set of terminals which can be connected to the terminals of other components, where the sets partition . Each connection is a pair of terminals, i.e. .
Iiia Connections
A terminal represents a point of connection for an electrical component. We use the quantities of real power, reactive power, voltage and voltage phase angle ( respectively) to model the flow of power into a component through a terminal. These are vectors in order to capture each time step in the horizon. For convenience we use a parent vector to represent all variables for a terminal , where . When two terminals are connected together, , we pose the following constraints:
The first two constraints ensure the power that leaves one terminal enters the terminal it connects to. The second two constraints ensure that the terminals have the same voltage and phase angle. This duplication of variables is necessary in order to decompose the problem for the ADMM algorithm. To avoid confusion, recall that connections and terminals are different from lines and buses (see Fig. 1).
We rewrite these constraints as for , where is the appropriate diagonal matrix. Further, we define the function as the LHS of this constraint for convenience: .
IiiB Components
At a high level, each component has a variable vector , an objective function , and a constraint function , where . For a component , the vector includes all terminal variables for that component: .
In the following sections we describe at a lower level the models for a variety of components. When necessary we use to index vectors by time, otherwise we imply standard vector operations. The index where is used to represent the value of the variable at the beginning of the current horizon, which we assume is known.
IiiB1 Bus
A bus has a variable number of terminals which depends on how many other components connect to it. For example, a bus might be connected to a generator, a load and 3 lines for a total of 5 terminals. Regardless of the number of terminals, the constraints take the form:
The first two constraints are an expression of Kirchhoff’s current law (KCL) in terms of power flows. The remaining constraints ensure that all terminal voltages and phase angles are the same.
IiiB2 Line
A line is a two terminal component which transports power from one location to another, typically from bus to bus. We model a line as having a constant conductance , susceptance and maximum apparent power . The AC power flow equations are derived from Ohm’s law, where :
(1)  
(2)  
(3) 
These constraints are identical for each time step, so we have left out the indexing by time to improve clarity. These equations are nonconvex, so they are often either approximated or relaxed, as we will discuss further in Section VII.
IiiB3 Generator
A generator is a single terminal component which produces real and reactive power. In our formulation the generator has a floating phase angle and voltage. A generator has lower and upper real and reactive power limits such that and , a ramping rate and a quadratic cost function for generation costs:
where is a diagonal matrix and .
IiiB4 Shiftable Load
A shiftable load is a single terminal component used to model electrical loads like dish washers and clothes dryers. These loads must start running between an earliest and a latest start time: . To model this we introduce binary variables for the horizon. A value of indicates that the component starts at the given time. A component runs for a duration of consecutive time steps, during which it consumes a load of .
A convex relaxation of this component can be obtained by relaxing the integrality requirement: .
IiiB5 Battery
A battery is a single terminal component with stored energy . The charge and discharge variables transfer energy to/from the battery:
where the constant gives the efficiency of the battery and the net power into the component is given by . The second constraint above prevents a greedy discharge of power at the end of the horizon.
Iv Alternating Direction Method of Multipliers
The goal is to minimise the sum of all component cost functions, subject to component and terminal connection constraints. This is a utilitarian view of the problem.
The augmented Lagrangian applied to the connection constraints is given by:
where is a penalty parameter and are the dual variables for the connection constraints. The dual variables also represent the locational marginal prices in our problem.
The alternating direction method of multipliers (ADMM) is a variation of the standard augmented Lagrangian method that enables problem decomposition [11, 14, 15]. A single iteration consists of two phases followed by a dual variable update. Components are each allocated to one of the two phases. The component sets and , and the variable vectors and represent this allocation.
The connections are split into three parts: , and . The intraphase connections () are those that are between components in (). The interphase connections are those where one component is in and the other is in .
The superscript is used to indicate the th iteration. At the start of the algorithm all terminal and dual variables are initialised to some values and . For the th iteration ADMM proceeds as follows:

Optimise for , holding constant at its value

Optimise for , holding constant at its value

Update the dual variables
For our optimisation problem this becomes:
(4)  
(5)  
(6) 
In the simple case when is constant, and are convex, and is affine, ADMM converges to a global optimum [11].
If a component has no intraphase connections, then its subproblem can be solved independently of all other components in the same phase. We adopt the partitioning scheme where contains all buses and the rest of the network. This allows us to fully separate all components, since buses will never connect to other buses () and nonbus components will never connect to other nonbus components (). As an additional benefit, some components are simple enough that we can analytically calculate the solution to their subproblem at each iteration, instead of invoking an optimisation routine. We adopt such an analytical approach for buses as proposed in [1]. Other partitioning schemes are discussed in Section VII.
Iva Residuals and Stopping Criteria
As in [1], we use primal and dual residuals to define the stopping criteria for our algorithm. The primal residuals represent the constraint violations at the current solution. We combine the residuals of all connections into a single vector . By indexing into the interphase connections , the primal residual is given by:
The dual residuals give the violation of the KKT stationarity constraint at the current solution. We collect the dual residuals for each each connection into the vector . For ADMM, the dual residuals are given by (see [11] for derivation):
These residuals approach zero as the algorithm converges to a KKT point. We consider that the algorithm has converged when the scaled norms of these residuals are smaller than a tolerance : , . Here is the total number of interphase terminal constraints minus the number of terminal constraints that are trivially satisfied (e.g., floating voltages and phase angles for generators). It is used to keep the tolerance independent of problem size.
V Implementation
We developed an experimental implementation of the ADMM method in C++ using Gurobi [16] and Ipopt [17, 18] as backend solvers. Gurobi is used for mixedinteger linear or quadratically constrained problems, and Ipopt for more general nonlinear problems. CasADi [19] was used as a modelling and automatic differentiation front end to Ipopt. This implementation was designed with flexibility in mind, so that a wide range of experiments could be conducted. For a more specialised and performance orientated implementation see [1]. Although this is a sequential implementation of the ADMM method, we timed the slowest component at each iteration to get an idea of how long a fully distributed implementation would take.
The experiments were run on machines with 2 AMD 6Core Opteron 4184, 2.8GHz, 3M L2/6M L3 Cache CPUs and 64GB of memory.
Vi Test Microgrid
The distributed algorithm we have presented can be used to control all levels of the electricity system, potentially in a hierarchical manner [1]. We focus our attention on a microgridlike distribution network since that is where we expect to have the greatest impact with the adoption of distributed generation and storage.
Our experiments are based around a modified 70 bus 11kV benchmark distribution network [20], which was chosen due to its comparable size to a suburb. We close all tie lines in the network in order to change it from a radial to a meshed configuration. This is because we expect microgrids to take on more of a meshed network structure to improve reliability, efficiency and to better utilise distributed generation.
The benchmark comes with a static PQ load at each bus, which we replace with a number of houses, depending on the size of the load. The houses are connected directly to the 11kV buses as we have no data on the low voltage part of the network. We assume that the power bounds we place on each household will be sufficient to prevent any capacity violation of the low voltage network.
A house is an independent agent that manages subcomponents. These include an uncontrollable background power draw, two shiftable loads, and optional battery and solar PV systems. A house has a single terminal through which it can exchange real and reactive power with the rest of the network. Each house has an apparent power limit of kVA.
We develop a typical house load profile by modifying an aggregate Autumn load profile for the ACT region in Australia (data from [21]). We assume that households consume on average 20kWh per day. This provides the basis for all uncontrollable household background loads. For the purposes of these experiments, we assume that the static PQ loads in the benchmark were recorded when loads were at 75% of their peak. We divide the benchmark static real power at each bus by how much power a typical house consumes at 75% of its peak power (1.45kW). Rounding down this number gives us an estimate of the number of houses which would be located at a given bus. This approach produces a total of houses for the network, about the size of an Australian suburb.
We place two generators in the network where the distribution system connects to upstream substations. These can be thought of as either dispatchable microgenerators or as representing the cost of importing power into the microgrid.
We randomise some of the generator and household load parameters to produce different problem instances, as can be seen in Table I. The time horizon spans 24 hours with 15 minute time steps, which produces a problem instance with over 2 million variables per horizon. The experiments were run with a primal and dual stopping tolerance of and a fixed penalty parameter of . To improve numerical stability, we scale the system to a perunit representation with base values at 11kV and 100kVA. This means that a real power residual of translates to 10W for a connection, or about 1% of the average household load.
Component  Param  Value  Units 

Generator  $/MWh  
Generator  $/MWhMW  
Generator  kW  
Generator  kvar  
House  kW  
House  kvar  
Shift 1  min  
Shift 1  kW  
Shift 2  min  
Shift 2  kW  
Shift 
In addition to the 70 bus network described here, we also ran a series of experiments on randomly generated networks similar to those described in [1]. These randomly generated networks ranged in size from 20 to 2000 buses, and were designed to be highly congested. We will occasionally mention some of the results from these random networks when they differ from our 70 bus microgrid.
Vii Impact of Power Flow Models
In this section we investigate how the ADMM method performs with different power flow models. We assess 4 alternative models to the one presented in [1], all of which give more accurate results, but are slower to converge.
Viia Power Flow Models
Due to their nonconvex nature, the AC power flow equations (1–3) are often either relaxed or approximated. Convex relaxations include a quadratic constraint (QC) model [22], a semidefinite program [23], the distflow (DF) relaxation [24, 25] and an equivalent SOCP relaxation [26, 27]. Approximations include the linear DC (DC) model [28, 29], the LPAC model [30] and the equations proposed by Kraning et al. (K) [1]. These alternative models are often used to solve difficult power network optimisation problems, for example, OPF, OPF with line switching, capacitor placement and expansion planing. ADMM was used with the DF model in [31] in order to control the reactive power of inverters on a radial network, and the K model was used with ADMM in [1] to solve the DOPF problem.
What is lacking is an understanding of the relative strengths and weaknesses of line models when used in a distributed algorithm. In this section we compare how the distributed ADMM algorithm performs when using the AC, QC, DF, DC, and K line models. We compare the differences in solution quality, feasibility, processing time and number of iterations for our test network.
ViiB Experiments
We generated 60 instances of our test network with the binary variables for the shiftable devices relaxed. To evaluate the quality of the solutions we calculate the percentage difference in objective value relative to the best known AC solution: . The means and standard deviations of the 60 instances are:
QC  0.031% (0.008%)  DF  0.039% (0.018%) 

DC  3.541% (0.072%)  K  4.726% (0.090%) 
We don’t have a guarantee that the AC solution is optimal because the equations are nonconvex. However, the relaxed models (QC, DF) provide a lower bound on the global optimum, allowing us to identify in the worst case how far the AC solution is from optimality. The DF model produces results that are slightly above the AC model, when it should be providing a lower bound. This is because the result is within the margin of error for our stopping criteria, which we estimated to be 1%. By decreasing we found the positive difference to shrink further into insignificance. Based on our results the optimality gap is tiny for both the QC and DF models. These results give us confidence that the nonconvex AC model, which is the only one that guarantees Ohm’s law is satisfied, produces solutions that are very close to optimal. Other work has come to a similar conclusion in a different setting [22, 13].
The DC model underestimates the optimal value by around 3.5% while the K model overestimates it by around 4.7%. Part of the reason for this is that the DC model ignores line losses and the K model is overestimating line losses.
Table II provides the number of iterations and time taken to converge in the form of means and standard deviations. The parallel solve time is the amount of time required to solve the problem in a fully distributed implementation. This was measured by summing together the time of the slowest component at each iteration. In absolute terms the parallel solve time is relatively small despite the fact that our implementation was designed with flexibility in mind, not performance. That said the K model is significantly faster relative to the other models. It converges in half the number of iterations required by the next fastest model, but as we have seen it gives us an inaccurate result.
Iterations (std.)  Time in sec (std.)  

AC  1945 (17)  148 (12) 
QC  1951 (14)  546 (33) 
DF  1933 (26)  110 (8) 
DC  4140 (50)  244 (8) 
K  1027 (52)  15 (1) 
The highly congested random networks that we generated produce similar results. However, for a number of instances the K model was infeasible (would not converge) where we had a valid AC solution, due to the exaggeration of line losses.
Fig. 2 shows an example of the primal residual convergence for different line models (the dual residuals are similar). The AC, QC and DF models overlap. One unintuitive result is the fact that the DC model converges poorly when it is in fact a very simple linear model. Large oscillations build up during the solution process which slows the rate of convergence. We performed a series of experiments in order to get a better understanding of this effect. The best explanation we have is that the DC model behaves like an undamped system, as it has no line losses and only linear constraints. The DC model will respond stronger for a given change in its terminal dual variables. The net effect is that oscillations build up across the network during the solving process. On the other hand the K model overestimates line losses, which means it is much less sensitive and no oscillations form. The AC, QC and DF models are somewhere in between these two extremes.
In reality communication delay will play a major part in the total solve time for the algorithm. For example, if we assume that it takes 60ms to communicate messages on our network (which happens a minimum of two times per iteration), then 1000 iterations would require 2 minutes of communication time. In certain circumstances it may be beneficial to cut down the total number of iterations, even if it requires more processing time per iteration.
We experimented with clustering some lines and buses together instead of fully decomposing the network. The idea is that, depending on the communication delay, this could reduce overall solve time by solving more of the problem at once. This was clearly beneficial for the random networks, but for the 70 bus microgrid the results so far are inconclusive. For example, clustering the 70 bus network into 16 parts halved the number of iterations required, but the overall solve time worsened. We expect this to change once we start to optimise our implementation, by reducing solver overhead and decomposing network components across time steps.
It is important to point out that we are giving the algorithm a naive starting point for both the primal and dual variables. In practice, the receding horizon control scheme will provide an excellent warm starting point, because the values from the previous horizon can be used for all but one time step. As a sanity check, we performed warm starting experiments for the AC model. Similar to what was done in [1], we duplicate a problem instance and then randomly resample the household background power and shiftable device power parameters according to the rule: . We used the solution of the original instance as a starting point for the modified instance. For the warm started run only needed 11% of the original iterations on average. In a second experiment we fully correlated the resampling step, which could represent a correlated change in solar panel output for many households. With , only 29% of the original iterations were required on average.
These results show the feasibility of using the nonconvex AC power flow equations for solving a distributed DOPF problem. The K model converges much faster, but it is unlikely to be usable in a realistic setting, as it ignores voltages and reactive power, and produces overly high costs. For these reasons, we adopt the AC model for the rest of our experiments.
Viii Discrete Decisions
We now extend the algorithm so that it can deal with discrete decisions. Components with discrete decisions are be common in realistic household models and they are required to model more complicated generators and network switching events. We develop a framework that can accommodate discrete variables otherwise we could risk network violations through unpredictable household behaviour.
Viiia Methods
We investigate three tractable methods for dealing with integer variables which have no global optimality guarantees. Just as we did for the AC equations, we will compare our result to a lower bound in order to get an understanding of the optimality gap. We categorise these methods as:

Relax and price (RP)

Relax and decide (RD)

Unrelaxed (UR)
The RP and RD approaches are broken up into two stages. The first stage, also known as the negotiation stage, is common to both methods. Here all integer variables are relaxed and the distributed algorithm is run until convergence. At this point the integer variables may take on fractional values, and this solution gives a lower bound on the optimal. In the second stage each component makes a local decision in order to force any fractional values to integers.
Recall from Section IIIB4 that shiftable devices have a binary variable for each time step, only one of which can take on the value 1 to indicate the starting time. In the second stage of the RP method, each house performs a local optimisation to determine how to enforce integer feasibility of . We designed a range of cost functions which penalise a component if it changes its terminal values from those negotiated. For a given cost function each house solves a MixedInteger Program (MIP) to obtain an integerfeasible solution. The two most effective cost functions that we identified are:
(7)  
(8) 
Where for a given house to bus connection is the negotiated terminal values for the bus and the negotiated dual variables. We use to represent the diagonal matrix where and is a penalty parameter.
The first function charges households at the negotiated price for what they actually consume, but they are also charged a quadratic penalty for operating away from the negotiated consumption. The second function requires the household to pay for all power that was negotiated in the first stage. Like the first function a penalty is charged for operating away from the negotiated operating point, however the penalty is scaled by the dual variables which can vary with time.
After this local optimisation step, we check that the solution is feasible and what the overall cost is. In order to do this we need to put some degrees of freedom back into the problem. In power networks the dispatch of generators are established in advance in response to an estimated demand. This forecast is never perfect, so a certain number of generators are paid to perform frequency regulation in order to balance demand in real time. In our experiments we employ both our generators for this use by allowing them to adjust their output. For these experiments we assume the same cost function and prices for both dispatch and frequency regulation.
In the second stage of the RD method, the largest value of each shiftable component is chosen to be fixed at 1 and the rest set to 0. After this is done the distributed algorithm is started again in order to converge to a new solution that is integer feasible.
The final approach, UR, attempts to enforce integrality satisfaction at each iteration of the distributed algorithm. We have already foregone theoretical convergence guarantees by our adoption of the nonconvex AC equations. Here we push the ADMM algorithm even further by allowing discrete variables into the algorithm (4–6), where Gurobi solves MIPs for houses, and Ipopt NLPs for lines.
ViiiB Experiments
We ran experiments on 60 random instances of our test network with a penalty of . The results can be seen in Fig. 3 which gives the fractional change with respect to the relaxed problem. This is shown both in terms of the cost of generation and the charge to households. The charge is the sum of household objective functions, which represents the amount of money they pay for their electricity. For the RP methods this is given by the cost functions in the previous section. For the RD and UR methods the charge is simply the final for each house.
For the relaxed problem itself, the true cost of generation can be different from the amount households are charged, as we are dealing with marginal prices. In addition to this, network congestion typically generates additional revenue above the cost of generation itself. An increase in cost for the integer feasible solution relative to the relaxed problem is an indication of the additional cost to the generators for balancing supply. Where household charge has increased relative to the relaxed solution, then households have decided to take on additional costs in order to change their consumption.
All methods produce costs that are within 1% of the relaxed problem. There is no significant difference between the four methods as they reside within our estimated margin of error based on our stopping tolerance. What these results suggest is that we have a tight relaxation of the integer problem. From other experiments, the relaxation of the shiftable devices is tight except where we have a heavily congested system. This combined with the fact that each shiftable load only contributes a tiny amount to the overall power demand gives us a tight relaxation.
The RP method only marginally increases solve time above the results in the Section VII. The RD method requires some of extra time as it performs a warm restart of the distributed algorithm. The UR method took 1.7 times longer on average as it solves MIPs during each iteration.
The charges to households are significantly higher for the RP method without gaining any benefit in terms of reduced costs. We ran the same experiments with a much smaller which all but eliminated charges, without any increase to costs. This suggests that for the sole purpose of managing discrete decisions, there is no need to have a strong penalty.
All of the methods we have presented provide an efficient means for dealing with the discrete decisions in a household. However, we will come to see in the next section why we favour the RP approach.
Ix Pricing Uncertainty
In this section we investigate the inclusion of stochastic components into our system. Many parameters such as background household power consumption, solar PV output can only be estimated. This means that the solutions that are negotiated through our distributed mechanism may not be applicable when it comes time to act. For example, if the output of household solar PV systems is lower than expected, then certain lines in the network could be overloaded if the network was already running near capacity.
Part of the way of dealing with this is to run receding horizon control so that we only ever act on the most immediate time step before reoptimising. However even within the most immediate time step we still have to deal with uncertainty. In this section we use the RP method from the previous section to encourage households to correct for any local sources of uncertainty. Of the three methods designed to handle discrete variables, this is the only one that can provide this type of incentive. To simplify the experiments we don’t perform full receding horizon control. With receding horizon control we expect the same trends, but with further improvement to costs.
Ixa Experiments
We perform a simple experiment on the 70 bus microgrid where we have added 2kW PV systems and 2kWh batteries independently and at random to half of the houses. The battery efficiencies were uniformly sampled from the interval . For normalised solar irradiance we use the simple relation: . We solve the first of the RP method with this irradiance, and then either lower or raise it it 20% before running the second stage. Fig. 4 shows the resulting costs and charges relative to the first stage result.
Lowering solar output produces increases to generation costs of between 810% relative to the original solution. Function 0 performs around 1% better on average and it increases the charge to the house in line with the increased generation costs. For function 3 the households are barely penalised for this increase in generation costs. When solar output is raised function 3 doesn’t react at all because it has no incentive to. Function 0 on the other hand takes advantage of the excess solar to lower overall costs by around 8%.
As we found in the previous section, generator costs are relatively independent of the RP cost function penalty parameter . However, can be used to increase household charges. This means that can be tuned over time to ensure that the market is budget balanced. For example, these extra penalties can be used to recover frequency regulation costs or to help pay for fixed network costs.
Larger penalties could also help to deter agents looking to game the system. In the current mechanism, there is nothing that prevents an agent from misreporting their behaviour during the negotiation phase. A larger penalty forces them to act on what they negotiated, so can significantly limit this behaviour. The problem is that it will be difficult to distinguish between someone looking to cheat the system, and someone with a large amount of local uncertainty or discrete loads. This makes it difficult to penalise one of these effects without penalising the others. Houses with batteries have more flexibility in how they can deal with local sources of uncertainty, so a system with batteries and reasonable penalties might provide the best compromise.
X Conclusion and Future Work
We have improved on existing distributed ADMM based DOPF methods by including more accurate line models and a twostage approach to manage discrete variables and uncertainty. We developed a suburb sized test microgrid, and found that the full nonconvex AC equations produce close to optimal solutions in short solve times. Our twostage approach provides a simple but effective means of managing household discrete variables and uncertainty in the network.
Future research will focus on investigating alternative distributed solving techniques with the aim of further improving convergence. There are also opportunities to parallelise the problem more by decomposing some components across time. Further research needs to be done on the realtime control that takes place within and between timesteps, and it might be possible to build a frequency regulation market into our distributed algorithm.
We need further experiments to investigate if our results carry over to larger discrete decisions, for example, those related to large industrial plant, generator startup costs, and line switching. We also plan to answer the important question of how susceptible this mechanism is to gaming in practice, and if this is a problem, what can be done about it.
References
 [1] Kraning, M., Chu, E., Lavaei, J., Boyd, S.: Dynamic network energy management via proximal message passing. Foundations and Trends in Optimization 1(2) (2014)
 [2] Gatsis, N., Giannakis, G.: Residential load control: Distributed scheduling and convergence with lost ami messages. Smart Grid, IEEE Transactions on 3(2) (June 2012) 770–786
 [3] MohsenianRad, A., Wong, V., Jatskevich, J., Schober, R., LeonGarcia, A.: Autonomous demandside management based on gametheoretic energy consumption scheduling for the future smart grid. Smart Grid, IEEE Transactions on 1(3) (dec. 2010) 320 –331
 [4] Ramchurn, S.D., Vytelingum, P., Rogers, A., Jennings, N.R.: Agentbased control for decentralised demand side management in the smart grid. In Tumer, Yolum, Sonenberg, Stone, eds.: Proc. of 10th Int. Conf. on Autonomous Agents and Multiagent Systems  Innovative Applications Track (AAMAS 2011), Taipei, Taiwan (May 2011) 330–331
 [5] MohsenianRad, A.H., LeonGarcia, A.: Optimal residential load control with price prediction in realtime electricity pricing environments. Smart Grid, IEEE Transactions on 1(2) (sept. 2010) 120 –133
 [6] Chen, Z., Wu, L., Fu, Y.: Realtime pricebased demand response management for residential appliances via stochastic optimization and robust optimization. Smart Grid, IEEE Transactions on 3(4) (dec. 2012) 1822 –1831
 [7] Scott, P., Thiébaux, S., van den Briel, M., Van Hentenryck, P.: Residential demand response under uncertainty. In: International Conference on Principles and Practice of Constraint Programming (CP), Uppsala Sweden (sep. 2013) 645–660
 [8] Tischer, H., Verbic, G.: Towards a smart home energy management system  a dynamic programming approach. In: Innovative Smart Grid Technologies Asia (ISGT), 2011 IEEE PES. (nov. 2011) 1 –7
 [9] van den Briel, M., Scott, P., Thiébaux, S.: Randomized load control: A simple distributed approach for scheduling smart appliances. to appear in IJCAI 2013 (2013)
 [10] Shinwari, M., Youssef, A., Hamouda, W.: A waterfilling based scheduling algorithm for the smart grid. Smart Grid, IEEE Transactions on 3(2) (June 2012) 710–719
 [11] Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J.: Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning 3(1) (Jan 2011) 1–122
 [12] Kim, B., Baldick, R.: A comparison of distributed optimal power flow algorithms. Power Systems, IEEE Transactions on 15(2) (2000) 599–604
 [13] Phan, D., Kalagnanam, J.: Some efficient optimization methods for solving the securityconstrained optimal power flow problem. Power Systems, IEEE Transactions on 29(2) (March 2014) 863–872
 [14] Douglas, Jim, J., Rachford, H. H., J.: On the numerical solution of heat conduction problems in two and three space variables. Transactions of the American Mathematical Society 82(2) (1956) pp. 421–439
 [15] Gabay, D., Mercier, B.: A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Computers & Mathematics with Applications 2(1) (1976) 17 – 40
 [16] Gurobi Optimization, Inc.: Gurobi optimizer reference manual (2014)
 [17] Wächter, A., Biegler, L.T.: On the implementation of an interiorpoint filter linesearch algorithm for largescale nonlinear programming. Mathematical Programming 106(1) (2006) 25–57
 [18] HSL Archive: A collection of fortran codes for large scale scientific computation (2014)
 [19] Andersson, J.: A GeneralPurpose Software Framework for Dynamic Optimization. PhD thesis, Arenberg Doctoral School, KU Leuven, Department of Electrical Engineering (ESAT/SCD) and Optimization in Engineering Center, Kasteelpark Arenberg 10, 3001Heverlee, Belgium (October 2013)
 [20] Das, D.: Reconfiguration of distribution system using fuzzy multiobjective approach. International Journal of Electrical Power & Energy Systems 28(5) (2006) 331 – 338
 [21] Australian Energy Market Operator (AEMO): www.aemo.com.au.
 [22] Hijazi, H.L., Coffrin, C., Van Hentenryck, P.: Convex quadratic relaxations for mixedinteger nonlinear programs in power systems. NICTA Technical Report http://www.optimizationonline.org/DB_HTML/2013/09/4057.html (March 2014)
 [23] Bai, X., Wei, H., Fujisawa, K., Wang, Y.: Semidefinite programming for optimal power flow problems. International Journal of Electrical Power & Energy Systems 30(6â7) (2008) 383 – 392
 [24] Farivar, M., Clarke, C.R., Low, S.H., Chandy, K.M.: Inverter var control for distribution systems with renewables. In: SmartGridComm, IEEE (2011) 457–462
 [25] Baran, M., Wu, F.: Optimal sizing of capacitors placed on a radial distribution system. Power Delivery, IEEE Transactions on 4(1) (1989) 735–743
 [26] Jabr, R.: Radial distribution load flow using conic programming. IEEE Transactions on Power Systems 21(3) (Aug 2006) 1458–1459
 [27] Bose, S., Low, S.H., Teeraratkul, T., Hassibi, B.: Equivalent relaxations of optimal power flow. CoRR abs/1401.1876 (2014)
 [28] Schweppe, F., Rom, D.: Power system staticstate estimation, part ii: Approximate model. power apparatus and systems, IEEE transactions on (1) (1970) 125–130
 [29] Stott, B., Jardim, J., Alsac, O.: Dc power flow revisited. Power Systems, IEEE Transactions on 24(3) (Aug 2009) 1290–1300
 [30] Coffrin, C., Hentenryck, P.V.: A linearprogramming approximation of ac power flows. Informs Journal on Computing (May 2014)
 [31] Šulc, P., Backhaus, S., Chertkov, M.: Optimal Distributed Control of Reactive Power via the Alternating Direction Method of Multipliers. ArXiv eprints (2014)