Surrogate models in bidirectional optimization of coupled microgrids

2.1 Upper-level model: Energy exchange

We assume that we have Ξ∈N many MGs which are partially coupled via transmission lines and can be interpreted as nodes of a non-complete graph, see Figure 1 where Ξ=4. Each MG κ, for κ∈[1:Ξ], consists of Iκ∈N agents modeled in detail in Subsection 2.2. We assume that each MG κ has an average power demand z¯κ(n) at time n. Given this, we can compute the total power demand Iκz¯κ of a MG. The control goal is to exchange power among the MGs in a way such that a desired quantity ζ¯(n) is targeted by controlling the residential storage units. We specify ζ¯(n) in Subsection 2.2, and assume for the moment that this quantity is known to the grid and has advantages for the grid operation. We assume that this desired quantity is independent of κ, but this is not necessary for the rest of the discussion.

Let δνκ describe the percentage of power Iνz¯ν(n) that is transferred from MG ν to MG κ. We enforce δνκ equals zero if there is no transmission line between the two MGs. Otherwise, the power demand of a MG κ is given by its own total power demand δκκIκz¯κ, where δκκ is what remains at the MG, and the sum over the power received from connected MGs, ∑ν≠κδνκIνz¯ν. For each time step in our prediction horizon of length N∈N≥2, we want to match this to the desired power demand starting at time k for N timesteps of each MG in a least-squares sense. The objective function is thus given by J:RΞN×RΞ×Ξ×N→R,

Here, the vector z¯=(z¯(k),…,z¯(k+N−1)⊤)⊤ with z¯(·)∈RΞ stacks the average power demand per MG and time step while the matrix η=(ηνκ)κ,ν=1Ξ∈[0,1]Ξ×Ξ incorporates efficiencies along the transmission lines.

We are interested in minimizing (1) under the following constraints: All exchange rates δνκ are within the interval [0,1], sum up to 1, and only transfer power in one direction, meaning that either δνκ or δκν is zero. Moreover, note that at this grid level, the average power demands per MG are known. Following [17], the optimization problem of the upper-level is, thus, formulated over the exchange rates δ,

The efficiency of a transmission line does not depend on the direction of the transfer, i. e., the matrix η is symmetric. Furthermore, we assume no loss without transport, i. e., ηκκ≡1 for all κ∈[1:Ξ], in the rest of the paper.

2.2 Lower-level model: Single microgrid

As we have seen in the previous section, we consider an average power demand at each MG as well as some desired quantity ζ¯(n). In order to understand these quantities better, we explain the modeling of the MG in all detail. The basis of our considerations forms the model presented in [24] and its extension [7]. Therefore, let each subsystem be equipped with an energy generation device (e. g., roof top photo-voltaic panels) and some storage device (e. g., a battery). Then, the i-th system, i∈[1:Iκ] in MG κ, can be described by the discrete time system dynamics,

Figure 2

Lower-level model: Star-shaped network of residential energy systems and central entity (CE). The quantity w=ℓ−g in (3b) is obtained from the data set [25].

The system can be controlled by charging uκi+ and discharging uκi− the battery at each time step. The length of a time step in hours is denoted by T>0, e. g., T=0.5 corresponds to a half-hour time window. The constants ακi,βκi,γκi∈(0,1] represent efficiencies w. r. t. self discharge and energy conversion. Furthermore, the initial SoC xκi(k)=xˆκi, where k∈N0 denotes the current time instance, is assumed to be known. State and input are subject to the inequality constraints,

For a concise notation we introduce the set Xκi=[0,Cκi] of feasible states, the set Uκi={(uκi+,uκi−)⊤∈R2|(4b)−(4d)hold} of feasible control pairs and the set

of feasible outputs over the next N time steps, i∈[1:Iκ], κ∈[1:Ξ]. Referring to the feasible sets of a MG κ we use the Cartesian product, e. g., D(κ)=Dκ1×…×DκIκ and z(κ)∈D(κ). Note that the sets Dκi, i∈[1:Iκ], and hence D(κ) are non-empty, compact, and convex.

The output quantity in (3b) is the power demand zκi of an individual agent in MG κ. The average power demand z¯κ=1Iκ∑i=1Iκzκi in each MG can then be computed from the individual power demands, and is used as an input to (1). We define ζ¯ in (1) as a stable reference trajectory by averaging over a past time horizon and over all individual residential units of all microgrids, the so-called overall average net consumption as proposed in [7],

where I=∑κ=1ΞIκ denotes the total number of agents within the entire smart grid. Due to averaging, the trajectory ζ¯ has little fluctuations and yields advantages for the grid operation.

Let us for the moment ignore the coupling described in Subsection 2.1. Then, δ(n) equals the identity for all n∈[k:k+N−1] and J becomes

Therefore, the overall objective J can be decoupled yielding the local optimization problem

per MG with local objective function g:RN→R≥0,

2.3 Fully coupled optimization problem

We are interested in optimizing the function (1). This function, in general, depends on δ as well as on z¯. As seen in the previous section, the average power demand z¯=z¯(u) depends on the control u, which we have to find in such a way that J is optimal. The overall optimization problem can be written as

Note that due to constraint (2c) the optimization of J w. r. t. δ is non-convex. Furthermore, the large scale of the optimization w. r. t. z¯ causes the use of a centralized solver to be expensive. In addition, using a centralized solver assumes the existence of a global entity gathering the information of the whole grid, in particular the personal data of each household, which is undesirable in practice. Hence, solving (7) centralized is impractical. In the subsequent section we present an approach to tackle (7) by solving the upper and lower-level problem iteratively. Doing so we avoid a node with full knowledge in the grid and only communicate specific aggregated information in each iteration among the agents.

3.1 Bidirectional optimization scheme

Assume that each MG κ, κ∈[1:Ξ], within the smart grid has already solved its inherent optimization problem (6) and based on the corresponding solutions z¯κ an energy exchange policy δ⋆ has been computed according to (2). This exchange yields an updated power demand

and hence, the difference Δz¯κ(n)=z¯κ(n)−z¯κ+(n) in power demand for all MGs. We are interested in updating z¯ in such a way that our cost function (1),

is minimized further. One could think about fixing δ and finding an optimal z¯. This, however, leads to an optimization problem not avoiding communication and coupling all microgrids. The trick here is now to fix not only the δ but also the z¯-components from all the microgrids but one. This leads to optimizing Iκ2(ζ¯+Δz¯κ)−z¯κ22 locally in each MG, where Δz¯κ is computed by using the z¯κ’s and δ’s from the previous optimization step. Intuitively, the difference Δz¯κ, κ∈[1:Ξ], can be interpreted as an additional load or generation, and, therefore, as a change of the desired power demand profile for MG κ. This yields the modified lower-level optimization problem

where ζ¯κ+=ζ¯+Δz¯κ In this formulation the updated reference trajectories ζ¯κ+, κ∈[1:Ξ], differ among the single MGs and depend on a given δ and a given previous z¯κ. The solution of the newly derived lower-level optimization problem can be solved with ADMM for all microgrids independently and in a parallel way.

Based on the updated reference value we solve (8) and (2) to improve the battery usage and the energy exchange and repeat the optimization until some terminal condition is satisfied, e. g., performance improvement less than a pre-defined tolerance or maximal number of iterations exceeded. This procedure is summarized in Algorithm 1. Note that we only update the reference ζ¯ on the lower level, since the upper-level optimization problem (1) does not change.

Algorithm 1

Iterative bidirectional optimization scheme

Neither convergence nor the interpretation of a potential limit of Algorithm 1 is clear a priori. Figure 3, however, experimentally shows convergence of the proposed scheme and a continuous improvement of the upper-level performance index. Here, we ran 10 iterations of the optimization scheme and plotted both the objective function values before and after the energy exchange within each iteration. The values stagnate after four iterations indicating that additional iterations do not further improve the overall performance. The next subsection elaborates on how to solve (8) in a fully distributed way using ADMM.

Figure 3

Evolution of the costs before and after energy exchange computed according to the bidirectional optimization scheme described in Algorithm 1, i. e., J(z¯j,δj−1) and J(z¯j,δj), resp. Note that J(z¯1,δ0) yields the costs without microgrid coupling.

3.2 Distributed optimization via ADMM

In this section we briefly discuss how to solve the lower-level optimization problem (6) or (8) using an Alternating Direction Method of Multipliers (ADMM) approach. We consider a single MG and therefore omit the index κ. Since the averaged output quantity appears in the objective function (6) or (8), we need to introduce an auxiliary variable a in order to decouple the lower-level optimization in the following way,

in a distributed way, cf. [5]. Following [7], the ADMM algorithm for (9) yields the three-step iteration ℓ↦ℓ+1,

According to Theorem 3.1 in [7] the optimization scheme (10) converges in the following sense.

According to [7] and the references [5, Section 3] and [4, Appendix C] therein, problem (6) fulfils the assumptions of Theorem 1.

4.1 Well-posedness

The following proposition states that for equality in (11), a proper mapping is defined. For a concise notation we replace the index κi by i here.

4.2 Radial basis functions approximation

Radial Basis Functions (RBFs) are used to interpolate functions based on a set of sampling data. We briefly recap some basics on RBFs. For a detailed introduction to theory and application see e. g., [20], for a similar approach where RBFs are used to replace ADMM we refer to [19].

Let M∈N denote the number of samples. Then, the interpolation function of (11) is given as the sum of basis functions ψm:RN×X×RN→R, m∈[1:M], and a regularization term q:RN×X×RN→RN. More precisely,

where χ=(w¯,x(k),ζ¯) is the joint inputs of Proposition 1, and αm∈RN, m∈[1:M]. The basis functions are so-called radial basis functions of the form, ψm(χ)=ψ(χ−χm), where the kernel ψ yields support close to the sampling data χm, m∈[1:M]. We choose an affine linear regularization q(χ)=β0+Bχ. Note that different choices are possible. The missing parameters αm, β0 and B are determined by interpolation conditions, cf. [19], [20].

In Figure 5, a possible fit via RBFs is visualized. Here, we interpolated given data from two-weeks of optimization (4540 data points) based on sampling data picking each 25-th data point to train (12). Then, we tested φRBF on the following day, and plotted the fitting. Our implementation is based on the Matlab toolbox DACE [27]. Note that the evaluation time of the RBF approximation grows with the number of data points used. Already with 180 data points to train (12) with N=6 causes the function evaluation of φRBF to be expensive, see Table 1. Using more data points would no longer yield an advantage over using ADMM w. r. t. computation time.

Figure 5

RBF and neural net fitting of the first component z¯(k) of z¯=(z¯(k),…,z¯(k+N−1))⊤ within 48 consecutive time steps indicating the quality of the approximation (11).

4.3 Neural networks approximation

Neural Networks (NNs) are a state-of-the-art method in artificial intelligence frameworks. Based on huge amounts of data M≫1 they are able to learn and recognize patterns in complex systems. We consider a NN of l-layers as an approximation to the mapping (11), i. e.,

where σ denotes the sigmoid function, and the weights W[l] and biases b[l] are determined during the training phase. Here, the number of neurons at layer l−1 and at layer l determine the number of rows and columns of W[l], respectively. Note that a separate neural network is trained for each MG. For an introduction to deep learning and neural networks we refer the reader to [28], [29]. To train the NNs, we used Matlab’s built-in toolbox nftool.

The overall goal of the approximation (13) is to be sufficient in the sense of the MPC performance shown in Figure 6. Our experiments in Figures 4 (bottom) and Figure 5 show that with one hidden layer of ten neurons only, a satisfying approximation on a 24-hours time window can be achieved if the training data is large enough. Note that NNs benefit from big data. In our case study, we trained the NN only on data corresponding to two weeks.

5.1 Model predictive control (MPC)

Consider the optimal control problem (6). In order to provide an optimal control sequence over an arbitrary long time horizon we use MPC. To this end, at current time instance k∈N0 we assume the future net consumption (wi(k),wi(k+1),…,wi(k+N−1))⊤∈RN to be predicted for all subsystems i∈[1:I]. Based on the prediction, Algorithm 1 is executed (per MG) to determine control sequences ui, i∈[1:I], and an exchange strategy δ. Then, only the first instances ui(k) and δ(k) are implemented and the time instance is incremented. Algorithm 2 outlines this MPC scheme.

Algorithm 2

MPC for coupled MGs

Note that Problems (6) or (8) and (2) have to be solved in order to determine z¯⋆ and δ⋆ in each MPC iteration. Therefore, the open-loop costs J(z¯⋆,δ⋆) can be computed in each iteration as well (cf. Figure 4). However, since Step 3 in Algorithm 2 suggests to only implement the first instance of the controls computed in Step 2, these costs are not attained. Instead the stage costs

are realized at each time step k∈N0 (cf. Figure 6).

5.2 Usage of surrogate models in MPC

We compare the performances using ADMM, RBFs, and NNs on the lower-level, i. e., in Step 4(b) of Algorithnm 1. In all numerical simulations we set T=0.5, N=6, Ξ=4, I1=50, and I2=I3=I4=10.^[1] The battery parameters were randomly chosen with mean values C=0.98, u_=−0.24, and u¯=0.25. Based on the battery capacities we set xˆi=0.5Ci. In order to incorporate losses along the transmission lines, we used the efficiency matrix,

For simplicity of the numerical computation, we only replaced the lower-level optimization routine for MG 1 and, thus, avoid training a separate surrogate model for each MG. We used Matlab for implementation.

Results on the MPC closed loop can be found in Figure 6 and Table 1. In Figure 6 the closed-loop performances of ADMM (black line) compared to perturbed ADMM, and ADMM (black line) compared to the two surrogate models are visualized. Similar to the open-loop case, small disturbances in ADMM have little impact and RBFs outperform the NN. The first column of Table 1 compares the sum of all MPC closed-loop performances using ADMM, RBFs and a NN while in column 2 the average runtimes of these approaches are reported. Note that when using a surrogate, we call ADMM once per MPC iteration. As elaborated in [7] in each ADMM iteration an N-dimensional vector has to be transmitted twice. Hence, both surrogates reduce the need for communication. Two great advantages of ADMM are that the local optimization (10a) can be parallelized and the global optimization is independent of the size of the MG. However, a single function evaluation such as (12) or (13) is faster than running the entire ADMM optimization routine.

Figure 6

Impact of mapping error (top) and approximation via RBF and NN (bottom) on the stage costs (14) within 48 consecutive time steps.

Note that in column 2 of Table 1 we ignored the communication between smart homes and CE which is needed to apply ADMM in practice. However, the runtime of ADMM impairs when executed in an actual smart grid while surrogates do not require additional communication.

In order to improve the performance of the NN, more sampling data has to be generated to increase the training set significantly. To avoid large offline computation times, we chose N=6, i. e., a prediction horizon of three hours, which is rather short compared to [7], [19].