Stochastic distributed model predictive control of microgrid with uncertain PV power prediction

This paper is concerned with a stochastic distributed Model Predictive Control (MPC) technique for power management of a photovoltaic (PV) generators-installed microgrid. The photovoltaic power supply has large uncertainty because it depends on weather conditions. To keep stable power supply to the microgrid, both accurate predictions of PV power supplies and efficient energy management based on the prediction are essential. We propose a distributed MPC method for microgrid management by combining the alternating direction method of multipliers-based distributed optimization and the randomized algorithm approach under the situation that a stochastic prediction model for the PV power prediction is available. The proposed method enables us to efficient energy management in a distributed way as well as the probabilistic guarantee of the line and battery capacity constraints. We demonstrate the effectiveness of the proposed method by a numerical simulation.


Introduction
A microgrid is a small-scale power network consisting of consumers, batteries and renewable energies. Recently, mass introduction of renewable energies has been expected from the environmental viewpoints. However, renewable generations have large uncertainty because they depend on weather conditions. To keep stable energy supplies to the microgrid, both accurate prediction of the renewable generations and efficient energy management based on the prediction are important. It should also be noted that renewable energies are introduced in various geographical locations, and it is difficult to manage them with the conventional centralized management. In the last decade, there have already been reported some works which focused on distributed microgrid managements including a distributed optimal power flow [1], distributed algorithms for managements of energy resources and loads [2], and a minimization problem of generation costs in the microgrid [3].
A model predictive control (MPC) is one of the effective techniques to cope with the uncertainties of renewable generations. For instance, Zheng et al. [4] addressed an energy management of the on-connected microgrid. Ravichandran et al. [5] and Pflaum et al. [6] considered the uncertainties of operation plans of electric vehicles and power consumption, and proposed stochastic MPC techniques for home energy managements by invoking randomized algorithms. Tomisawa et al. [7] reported a stochastic MPC method of power distribution networks taking account of the temperature variations of distribution lines due to weather conditions. The references [8,9] addressed the stochastic MPC techniques for a simplified microgrid model in which consumers, PV generators, and batteries are consolidated into single agents. Especially, with the aid of the randomized algorithm [10], Takeda and Takaba [9] devised a stochastic MPC method of a microgrid with PV Power prediction based on the Just-In-Time (JIT) modelling technique [11].
For a more realistic microgrid model consisting of a number of batteries, consumers, and PV generators as well as battery management systems (BMSs), Takeda and Takaba [12] tried to extend the stochastic MPC technique proposed in [9] and applied the local MPC for each battery. However, this technique does not optimize the global benefit or performance of the whole microgrid. In the deterministic settings under the assumption that the PV generations are exactly predicted, Namba et al. [13] proposed a dual decomposition-based distributed MPC technique which guarantees the global optimal for the whole microgrid. Hassan et al. [14] proposed a distributed algorithm for the optimal power flow problem with the chance constraints for the PV power prediction by invoking the Alternating Direction Method of Multipliers (ADMM).
In this paper, we propose a microgrid management technique by manipulating charge-discharge behaviour of each battery in a distributed and cooperative fashion based on the combination of the ADMM and a randomized algorithm [10] under the situation that the stochastic PV power prediction [15] is available. One of the major contributions of this work is that the proposed method enables us to efficiently manage the microgrid in a distributed way as well as the probabilistic guarantee of the line and battery capacity constraints. To the best of author's knowledge, a similar stochastic distributed optimization technique for power systems was only proposed in the reference [14]. Hassan et al. [14] only addressed the static optimal power flow problem whereas we consider the dynamical energy management using the model predictive control. The dynamical control of the microgrid based on the receding horizon strategy is more practical because the PV power is heavily time-varying due to the weather variations as already mentioned in this section.
Preliminary results of this work were partially presented in the conference paper [16]. In this paper, the problem setting is improved to a more realistic one which takes account of electricity prices more appropriately. More numerical simulations are included to verify the effectiveness of our proposed method.

Solar irradiance and PV power prediction
In this section, we briefly review the JIT-based solar irradiance proposed by Suzuki et al. [11] and its modification [15]. It is known that the electric power generated by a PV panel is determined by the specification of the PV panel and the solar irradiance at the point of interest. Hence, we first discuss the solar irradiance prediction method based on the JIT modelling technique, and then review how convert the solar irradiance to the PV power. Interested readers are encouraged to refer to [11,15] for the JIT modelling approach to solar irradiance predictions and its modification.

Solar irradiance prediction via just-in-time modelling approach
The output prediction procedure of an unknown system for a given query input based on the JIT modelling is as follows. Firstly, we construct a database consisting of past input-output datasets. From the database, we select a small number of datasets whose input data are the nearest to the query input. These datasets are referred to as "neighboring datasets." The output for the query input is predicted from the neighbouring datasets based on the idea that outputs for close inputs are close to each other.
In the solar irradiance prediction, we employ weather forecast data and solar irradiance measurements as the input and output data, respectively. We obtain the weather forecast data from the Meso-Scale Model Grid Point Values (MSM-GPV) data of [17], and the solar irradiance measurements from [18].
We employ 11 meteorological factors for the solar irradiance prediction: atmospheric pressure, wind (south-north/east-west), temperature, humidity, rainfall, cloudiness (high/middle/low), solar attitude, and day length. Since the MSM-GPV data are updated every 3 hours, i.e. 8 times a day, we prepare 8 databases corresponding to the updating times. In the database at T o' clock, the meteorological factors are composed of hourly sequences from T o' clock to 23 o' clock. We execute the following steps at T o' clock, T = 0, 3, 6, . . . , 21.
(1) Obtain the weather forecast data of the day as the query input. (2) Select from the database K neighbouring datasets whose input data are the nearest to the above query input. (3) Estimate the future solar irradiance sequence by weighted averaging the past solar irradiances of the selected K datasets where y j T is the jth measured solar irradiance on the horizontal surface in the database at T o' clock, and j 1 , j 2 , . . . , j K are the indices of the neighbouring datasets, and a j κ ,T is the weight.
It may be noted that we update the JIT-based predictions every 3 hours according to the update of the MSM-GPV data, while the original method due to [11] is the one-day-ahead prediction carried out only once a day.

Modification of the JIT-based prediction
Since the conventional JIT-based irradiance prediction utilizes only the weather forecast data, its prediction accuracy may be degraded by weather variations during the day. On the other hand, the accurate on-line prediction of PV generations is required for the model predictive control of a microgrid.
For this purpose, we add an on-line correction to the JIT-based prediction using the most recent past measurements of solar irradiance.
Let y be the true solar irradiance, and define the prediction error at time t Note that T denotes the most recent updating time of the JIT-based prediction available at time t. Since y T has time correlation, we assume that y T is modelled by the autoregressive (AR) process model where a l,T , l = 1, 2, . . . , L are the regression coefficients, and w T is the Gaussian white noise process with mean 0 and the variance σ 2 T . The regression coefficients a 1,T , . . . , a L,T and the variance σ 2 T are identified by the least squares method using the neighbourhood datasets y j κ T , κ = 1, . . . , K. For example, in the case of the one-day ahead prediction at 0:00, the coefficients (a 1,0 , . . . , a L,0 ) are identified by where K AR is the number of the nearest neighbourhood datasets used for the parameter estimation. To achieve good prediction accuracy, we construct the AR model from only K AR neighbourhood datasets. Letŷ T (t + k | t) denote the k-step prediction of y based on the measurement up to time t. We propose the corrected prediction model of y T (t + k) aŝ where and the initial valuesŷ . . , L. Note that, we use the on-line corrected PV power prediction model (2) during the daytime because y(t) = y JIT (t) = y(t) = 0 obviously holds from sunset to sunrise.

PV power prediction
To obtain the PV power prediction, it is necessary to convert the solar irradiance on the horizontal surface to PV power generation taking account of the characteristics of PV panels such as conversion efficiency, tilt angles, etc. We approximate the relationship between the horizontal solar irradiance and the PV power generation by the following equation [19] and obtain the PV power prediction bŷ where S andŜ T denote the PV power and its prediction, W nmo is the nominal maximal output, PR is the performance ratio of the PV panel, and α is the conversion coefficient from the solar irradiance on the horizontal surface to that on the tilted PV panel surface.

Microgrid model
In this paper, we consider a microgrid model consisting of M consumers, M PV generators and M batteries.
Considering realistic situations, we add the following assumption.

Assumption 3.1:
(1) The supply-demand relationship between consumers and PV generators is one-to-one. (2) Consumers, PV generators, and batteries are not connected to any other consumers, PV generators, and batteries, respectively.
Notice that the above assumption describes the connections between the nodes by the distribution lines. This assumption guarantees that all the power flows in the microgrid model can be uniquely determined without any contradictions. An example in the case of M = 3 is shown in Figure 1 [13]. In this example, solid arrows show power flows, dashed arrows show information flows. F XY denotes the power flows from the node X to the node Y. A double-headed arrow means that bidirectional power interchange is allowed. Especially, F UD i denotes the power flow from the main grid to the ith consumer. In this paper, we separate the power flow from the main grid to the ith consumer as where F sell UD i ≤ 0 and F buy UD i ≥ 0 denote the selling and buying power flows of the ith consumer to/from the main grid, respectively. Each battery is equipped with a battery management system (BMS). Each BMS attempts to optimize a local benefit or objectives of its managed consumer (see Section 3.2). For this purpose, the BMS determines the power flows from and into the corresponding battery and the power flows between the main grid and the corresponding consumer. For example, the decision variables of the BMS1 are F B 1 D 1 , F S 1 B 1 , F sell UD 1 , and F buy UD 1 . We also assume that each BMS exchanges information with the coordinator that aggregates the information of the whole microgrid.
Under Assumption 3.1, there are various constraints in the microgrid as follows.

Battery capacity constraints
Line capacity constraints Charge-discharge dynamics of batteries where t is the sampling period, B i (t) is the charge of the ith battery, B max i is the maximum charge capacity of the ith battery, D i is the index set of the consumers connected to the ith battery, S i is the index set of the PV generator's connected to the ith battery, F XY is the minimum line capacity of the flow F XY , and F XY is the maximum line capacity of the flow F XY . In the case of F XY < 0, the reverse power flow from Y to X is allowed. We also assume that the ith PV power is delivered to the connected consumer and batteries, namely, where S i (t) denotes the ith PV power and B i denotes the index set of the batteries connected to the ith PV generator. In this paper, we employ the model predictive control method for the microgrid management. In the MPC strategy, we solve a certain finite horizon optimal control problem at each time instant, and apply the initial value of the optimal control sequence with shifting the horizon. To simplify the notations, we denote the k-step future prediction x(t + k|t) based on the data up to the current time t by x[k]. In particular, we have x[0] = x(t). We also abuse this notation to express x[k] := x(t + k) when the future value x(t + k) is exactly known a priori. Remark 3.1: From a viewpoint of the hierarchical structure of power systems, this paper addresses the high-level control of energy interchange/trading between the main grid and the consumers through charge-discharge manipulations of batteries. It is thus assumed throughout this paper that the lower-level control of voltage regulation has already been achieved by power electronic devices. It may also be noted that the high-and low-level controls are operated in the quite different time scales.

Formulation of distributed MPC for microgrid
We take the following control objectives for the purpose of stable energy supplies to the microgrid.
(2) Energy reservation of battery charges to cope with unanticipated faults. (3) Low purchasing costs and high selling profits through the power interaction between the consumers and the main grid.
The control objective 1, that is the primary objective of the microgrid management, is expressed as where D j [k] := D j (t + k) denotes the future demands of the jth consumer. For simplicity of discussion, we assume that D j [k] is exactly known a priori throughout this paper. It should be noted that our method can be easily extended to a more practical situation where D j [k]'s are uncertain and subject to prediction errors, if an appropriate stochastic prediction model for the demands is available.
To achieve the control objectives 2 and 3, we define the local objective function for the ith BMS by where N is the horizon length, and the positive constants P it , γ , Q sell , Q buy , and R i are the weights. In particular, Q sell and Q buy are the selling and buying unit prices of electricity of the consumers to the main grid, respectively. The first term of J i corresponds to the control objectives 2. The second term represents the control objective 3. The third term aims at smooth evolutions of power flows. In this paper, we choose the weights P it as where β i0 , β i1 , P 0 it are positive constants. It may be noted that P it is affine in the difference of the predicted PV powerŜ i (t) and the demand D i (t) of the corresponding PV generator and consumer, respectively. We can expect that the above weight enables us to efficiently manage the battery charging. For example, if the PV power is sufficiently larger than the demand, then the BMS behaves to store more energy from the PV generator. On the other hand, if the PV power is smaller than the demand, then the BMS discharges to the corresponding consumer to keep low dependency on the main grid under the supply-demand balancing.
Taking account of the practical situation in Japan [20], we make the following assumption.

Assumption 3.2:
The buying electricity price is always greater than the selling electricity price, namely, Q buy > Q sell .
Note that we can eliminate either one of F sell UD i and F buy UD i by using (7). For example, we can replace the second term on the RHS of (15) by Hereafter, we will only consider F buy UD i as a decision variable of the ith BMS.
Let us denote the decision variables of the ith BMS by x i , and the collection of all the decision variables of the BMSs by For example, the decision variable for the BMS 1 in Figure 1 is . . .
In this paper, we assume that the solar irradiance y(t) and PV power S(t) are expressed by (2) and (5), respectively. Correspondingly, we execute the solar irradiance prediction and the PV power prediction using (4) and (6), respectively. By the PV power prediction technique mentioned in the above, the objective function J i and the constraints expressed in (9), (10), (13), and (14) contain the random variable w for the predicted PV power. Let us introduce the following assumption. It is easily noticed that geographically close PV panels generate similar power outputs. By Assumption 3.3, we claim that such correlation among PV powers is represented by the JIT-based predictions, and the random variables w 1 , . . . , w M are the residuals which are uncorrelated with each other.
Let us introduce the chance constraint where ∈]0, 1[ is determined by the designer, and 1 − is the lower bound for the probability V(x) for some x. Let the inequality g(x, w) ≤ 0 represents the inequalities (9) and (10) which contain the random variables with the redundant variables eliminated by using (13), (14). We also denote the inequality constraints (8), (11), and (12) which do not contain the random variables by where X i is a convex and closed constraint set. We wish to solve the following optimal control problem (Chance-constrained Convex Optimization Problem) at each sampling time, in order to find the optimal power flows from and into the batteries, It may be noted that the global objective function for the whole microgrid is defined by summing up the local objective functions of BMSs.
Recall that, in the receding horizon strategy, we apply only the initial values (i.e. the values at k = 0) of the optimal inputs. In the example of Figure 1,

Stochastic distributed MPC algorithm
In this sub-section, we will derive a distributed algorithm for solving (CCOP) based on the randomized method [10] and the ADMM method [21,22]. For completeness of this paper, the basics on the convex optimization with the ADMM are summarized in Appendix.
Consider the objective function J i for the ith BMS defined in (15). The second term of the objective function may contain the other BMS's decision variables x j . On the other hand, from (13) and (14), the second term is linear with respect to x ι , ι = 1, . . . , M, and the random variable w i . Thus, we can rewrite J i as and this function is separable. It should be noted that J ii is a quadratic function of (x i , w i ), and J ij is a linear function of x j . Taking account of this fact and Assumption 3.3, we introduce a new objective function for the ith BMS as Thus we can express J, J i by the following quadratic functions The definition of J i ensures the positive semidefiniteness of the matrix H i , namely the function J(x) is not strictly convex but convex. Next, we consider the chance constraints in (CCOP). Since it is difficult to directly handle the chance constraints, we employ the randomized algorithm due to Calafiore and Campi [10]. In this paper, we generate ν random samples w (q) , q = 1, . . . , ν according to its probabilistic distribution which we got through the parameter identification for (2), and replace the chance constraints by the deterministic inequalities g(x, w (q) ) ≤ 0, q = 1, . . . , ν. Recall that g(x, w) is affine with respect to x and w. Then, we can express g(x, w) by where A i is a constant matrix and b i is a linear function of w i . Based on the above discussion, we relax the chanceconstrained convex optimization problem (CCOP) into the following sampled convex quadratic programming (SCQP1), The following fact holds true for (SCQP1) due to Corollary 1 of the reference [10].
is satisfied, where T is the number of the decision variables. Then, the optimizer x * for the sampled convex quadratic programming (SCQP1), if exists, satisfies the probabilistic inequality

Remark 3.2:
In this paper, we assume that the ith BMS only carries out the PV power prediction of the corresponding ith PV generator, and cannot access the power predictions of the other PV generators. We also assume that the ith BMS generates the random samples w (q) i , q = 1, . . . , ν for the random variables w i , i = 1, . . . , M contained in the ith PV power. The number of the samples is determined by the coordinator based on (16), and sent to all the BMSs. Recall that w i and w j , j = i are independent of each other by Assumption 3.3. It is possible that the coordinator carries out the predictions of all the PV generator's and the generations of all random samples. However, this imposes heavy computational duty on the coordinator. To avoid this difficulty, we do not employ it.

Remark 3.3:
The number of the inequality constraints in (SCQP1) increases according to the number ν of the random samples. To avoid this computational complexity, we extract the tightest constraint from ν inequalities, where b i is the vector consisting of the minimal components among b i (w (1) In the remainder of this section, we shall apply the ADMM method [21,22] to (SCQP1) and propose a distributed MPC strategy for the microgrid management. Note that the merit of the ADMM method is that the convergence speed is relatively fast and applicability to the optimization problem with the nonstrictly convex objective function. As widely known, the ADMM enables us to solve a certain class of convex optimization problems in a distribute way [21]. To reduce (SCQP1) to the form suitable for the distributed optimization, we introduce the non-negative slack variables ξ i and rewrite the inequality constraints as the equality constraints where A i is the coefficient matrix corresponding to x i and the the decision variables are redefined by x : We also introduce the box constraints on the slack variables ξ i as (18) where ξ i is a sufficiently large positive constant vector representing the upper bound on ξ i . The constraint (18) guarantees the satisfaction of Assumption A.1 (iii), and enhances the numerical stability of the proposed algorithm.
Let us now reformulate the deterministic sampled convex programing (SCQP1) as the equality constrained quadratic convex programming (SCQP2) which is in the separable ADMM form (SCQP2) minimize where z i , i = 1, . . . , M is the auxiliary variable. By invoking the techniques from Section 5.4.2 in the reference [21], the distributed ADMM algorithm for (SCQP2) consists of the following three equations.
where ω and λ are viewed as the incentive terms for the global optimality based on the local decision variables at the current iteration.

Algorithm 1 Stochastic distributed MPC
Step 0 The coordinator determines the step size c and the number of the samples ν based on (16), and sends them to BMSs.
Step 1 At the current time t, each BMS generates the ν samples of the random variable w i and extractsb i which contains the information of the tightest inequality based on (17), and send it to the coordinator. After the preparation steps described above, the distributed ADMM algorithm is carried out.
Step 1b At the current iteration τ , for given λ τ and ω τ , the BMSs solve the sub-problems (19) and determine the power flows from and into the batteries. Each BMS sends the x i,τ +1 to the coordinator.
Step 1d If the termination criterion is satisfied, then take x i,τ +1 as the global optimizer and proceed to Step 2. Otherwise, take τ := τ + 1 and go back to Step 1b. For the termination criterion, we employ ||ω τ || ≤ η, where η is a design parameter.
Step 2 Apply only the initial value of the predicted optimal control sequence to the BMSs. Then, shift the current time step to t := t + 1, and go back to Step 1.
In the framework of the microgrid management, the computation of (19) is carried out by each BMS, while the update of the penalties (20) and (21) are carried out by the coordinator. This means that the coordinator is monitoring and imposing the penalties for constraint violations on all the BMSs. Note that the iteration counter τ is different from the sample step k in the MPC strategy. During an intersampling interval, the BMSs and the coordinator exchange information to calculate (19), (20) and (21) iteratively. The stochastic ADMM-based distributed microgrid management is summarized in Algorithm 1, which solves the optimal control problem (SCQP2) at each MPC step. As stated in Fact A.1, the subsequence {x i,t } generated by the (19)-(21) asymptotically goes to one of the primal optimizers. It should also be noted that the proposed iterative algorithm needs to converge within the sampling period. The information structure of Algorithm 1 is shown in Figure 2.

Numerical simulations
In this section, we demonstrate the effectiveness of the proposed method through numerical simulations. To illustrate the effectiveness of the stochastic MPC, we carry out the deterministic MPC in which we take Though Hassan et al. [14] studied a chanceconstrained distributed algorithm for the optimal power flow problem, their method is not directly applicable to our problem, and we will not conduct a numerical simulation for comparison with [14]. This is because the reference [14] addresses the manipulation of active and reactive powers of inverters to achieve the optimal static power flow, whereas this paper is concerned with the charge-discharge manipulation of energy storage devices.

Simulation settings
We show the information of the PC and software used in the simulation in Table 1.
We employ the distributed microgrid model in Figure 1 for the simulation. We define the physical parameters as in Table 2 with reference to the small-scale simplified microgrid model due to Math-Works [23]. Taking account of practical situations, we set Q sell = 15 [JPY/kWh] and Q buy = 20 [JPY/kWh], respectively, due to the electricity price plans of Kansai Electric Power Company [20].
The profiles of the true PV generations and demands are depicted in Figure 3. The profiles of the weights P it are shown in Figure 4. In this simulation, we employ the on-line corrected PV power prediction model (2) from 5:00 to 19:00.

Simulation results
The results of our proposed method (stochastic MPC) are shown in Figure 5. Figure 5(a) shows the chargedischarge behaviour of the batteries and Figure 5(b-d) show the power flows related to each BMS based on the stochastic MPC technique. Dashed lines show the battery capacities and line capacities for the power flows shown in corresponding colours. As can be seen from Figure 5(a), from the 0 o' clock to the morning, we can see that all the batteries discharge to their corresponding consumers to keep the storages at certain levels. Moreover, the transmission of PV-generated power through the corresponding consumer makes each battery near full charge. From the noon to 15 o' clock, excessive PV power is sold to the main grid through the corresponding consumer. It should be noted that the power flows satisfy the constraints for at all sampling times. In the case of the deterministic MPC shown in Figure 6, although the charge-discharge behaviours of the batteries are similar to the stochastic MPC's, the power flows F UD 1 , F UD 2 and F UD 3 violate the line capacities at some sampling times because of the insufficient consideration of the PV prediction errors. We also carry out simulations for other 11 days with 10 trials for each so that verify the effectiveness of our proposed method. Notice that our proposed method does not guarantee the constraints satisfaction deterministically. The results are summarized in Table 3. In Table 3, "infeasible" means that the optimal control problem (SCQP2) became infeasible. The stochastic MPC-based approach satisfied the constraints for 11 days with all the trials except for one infeasible case, whereas the deterministic method satisfies only two days.
Next, we verify the economic effectiveness of our proposed method in comparison with the previous work [16]. Recall that the net electricity charge for the ith consumer is given by In [16], the objective function was used, where Q i is a positive weight. In the simulation, we choose Q i = 6, i = 1, . . . , M, with P it and R i the same as in the previous paragraphs. It should be noted that this objective function does not consider the trading loss/profit of electricity. To make a fair comparison, we set F buy UD i and F sell UD i as follows for the previous method of [16]. It can be shown that, for the power flow F UD i obtained the previous method, this choice of F sell UD i and F buy UD i minimizes the net electricity charge of (22).
The net electricity charges for the consumers and whole grid are summarized in Table 4. We can observe that the economic performance is improved on the extension in this paper. From the above discussions, it is seen that our proposed method is effective for microgrid management.

Concluding remarks
In this paper, we proposed the stochastic MPC technique for the distributed microgrid management by invoking the ADMM and the randomized algorithm [10] under the situation that the stochastic PV power prediction described in the Section 2 is available. The proposed method probabilistically guarantees the satisfaction of the constraints in the microgrid, such as the line capacities and battery capacities. The numerical simulation illustrates the effectiveness of our proposed method. Especially, we have shown that our proposed method can deliver the economic benefits to consumers.

Appendix. Convex optimization via ADMM
In this section, we give a brief review of the ADMMbased convex optimization due to Bertsekas [21] and Boyd et al. [22]. The ADMM is one of the iterative algorithms to solve convex optimization problems based on the augmented Lagrangian. It is known that the convergence speed is relatively fast for many practical examples [22].